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ABSTRACT 



Positron Annihilation Spectroscopy (PAS) has become a useful tool for 
the nondestructive Investigation of the defect structure of materials. Of the 
existing PAS measurement techniques Ooppler broadening appears to offer 
several practical advantages. These advantages Include low radioactive source 
strength requirements, comparatively short experimental times, and minimal 
geometrical restrictions. Unfortunately, the technique suffers from some 
Inherent disadvantages. These Include comparatively poor instrumental energy 
resolution and a limited ability to distinguish between multiple defect types. 

In an attempt to overcome the difficulties associated with the use of 
Doppler broadened PAS the Fast Fourier Transform/Power Spectrum (FFT/PS) de- 
convolution alyorixhm has been developed in this work. This algorithm is com- 
putationally efficient, nonlteratlve, and requires minimal a priori informa- 
tion. The utility of the FFT/PS deconvolution algorithm is demonstrated for 
high purity aluminum samples containing known concentrations of vacancies or 
prismatic dislocation loops. It is shown that the FFT/PS procedure 1s capable 
of automatically correcting for accidental and/or intentional changes in the 
instrument resolution function, In addition, based or, the limited data anal- 
yzed in this work, it appears that both the sensitivity and the accuracy of 
the energy distribution llneshape parameters are improved by deconvolution. 
Finally, it is demonstrated that if the filtering error introduced by the 
FFT/PS procedure is controlled, then 1t is possible to distinguish between 
samples containing predominately vacancy type defects and those with prismatic 
dislocation loops as the dominant defect type. 
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CHAf TER I 
INTRODUCTION 



Over SO years ago. Dirac postulated the existence of an atomic particle 
with a mass equal to that of an electron, but which is positively charged. In 
1932 these positively charged electrons, referred to as positrons, were first 
observed by C. D. Anderson as a by-product of cosmic ray showers. Since that 
time positron research has been pursued In such varied disciplines as high en- 
ergy nuclear physics, astrophysics, and materials science. It is within the 
broad context of materials science, and more specifically in the area of non- 
destructive testing (NOT) of aluminum, that the Interaction of positrons with 
matter will be examined in this work. Before going on to describe the use of 
positrons as an NOT tool in the investigation of the defect structure of 
metals, a brief review of positron physics will be presented. A comprehensive 
treatment of this subject has been written by West (1974). 

As postulated by Dirac In 1935, positrons and electrons are complemen- 
tary states of matter such that when the two particles encounter one another 
they will annihilate each other. Their kinetic and rest mass energies will be 
dissipated as the energy of one, two, or three gamma rays 1n such a way as to 
conserve momentum. It can be shown that the ratio of the two-photon to three- 
photon annihilation event Is 372:1 (Stevens 1980), while the ratio of the two- 
photon to one-photon events is ~2.6xlo6 : l, (West 1974). With these ratios as 
justification the remainder of this work will consider only the annihilation 
event yielding two photons. If both the electron and positron were at rest, a 
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positrons in vacancy type lattice defects. This effect was first observed by 
MacKenzie et al_. and Dekhtyar et aU In 1964. Schultz et a]_. (1981) describe 
the changes in the annihilation characteristics due to the localization of a 
positron In a raonovacancy type defect as follows: 1) Because the localized 
positron has a smaller overlap of its wave function with that of the high ener- 
gy core electrons than 1t does when diffusing through the nearly perfect lat- 
tice, the annihilation rate with high energy core electrons decreases. 2) The 
density of conduction electrons In the vacancy is also decreased resulting in a 
decreased annihilation rate with conduction electrons. However, the percentage 
reduction is less than in the case of core electrons so that the net result is 
an inc-^ase in the fraction of annihilation events involving low energy conduc- 
tion electrons, thus resulting in a decrease in the width of the Doppler broad- 
ened annihilation radiation distribution. 3) The potential for electrons in 
the vacancy is higher than normal resulting in the conduction electrons travel- 
ling more slowly. This in turn also decreases the observed Doppler broadening. 
4) A second order effect is that the overlap of the positron wave function with 
weakly bound outer-core electrons 1s reduced less than is the overlap with 
inner-core electrons. This too leads to less broadening of the observed 511 
KeV line. S) The uncertainty principle requires that localization of the posi- 
tron be accompanied by an increased spread in its momentum distribution (and a 
corresponding decrease in its de Broglie wavelength). As a result, the contri- 
bution of the positron to the center-of-mass velocity increases corresponding 
to a broadening of the annihilation energy spectrum (Schultz et aK 1981). The 
first four effects lead to a narrowing of the annihilation energy spectrum 
while the fifth effect tends to broaden the spectrum. The combination of these 
effects results in more sharply peaked (narrower) Doppler spectra for metals 
c .°.?l < - 1 n1 "« ^5fe cjts . whjch are . capable Qf_tra^p_in3 jjc^ ^ojisJ5j.e geJ__1 9ftn) . 
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As discussed by Siegel (1980) It Is Important to note that the posi- 
tron Is not equally sensitive to all lattice defects. It is sensitive to 
defect regions which can be characterized as having a lower-than-average elec- 
tron density, for example, vacancies, small vacancy clusters, voids, and dis- 
location lines or jogs. On the other hand, positrons are relatively insensi- 
tive to defect regions with a hlgher-than-average electron density such as 
interstltials. Indeed this defect specificity 1s one of the advantages PAS 
measurements offer as compared with other NDT techniques such as residual 
electron resistivity (Siegel 1980). 

In addition, positrons are capable of probing the entire depth of the 
valence band (MacKenzie et aU 1964) which is in contrast with techniques such 
as electrical and thermal conductivity, thermoelectric power. Hall coef- 
ficient, magnetic susceptibility, and Knight shift, all of which yield infor- 
mation restricted to the electrons near the Fermi surface. 

PAS measurements are an attractive NDT tool since the positron causes 
negligible perturbation of the electronic structure as it passes through a 
crystal (Johnson et al.. 1978) and the positron -electron Interaction does not 
appreciably alter the momentum distribution of the electrons (Szotek and 
Chmlelewskl 1981). A further advantage of the technique is that 1t yields 
information regarding both the electron momentum distribution as well as the 
electron density within the sample (Siegel 1980). 

The three most common PAS techniques are lifetime measurements, angu- 
lar correlation, and Doppler broadening. The effective electron density in 
the region from which the positron annihilates can most easily be determined 
using lifetime measurement techniques. These measurements usually use 22 Na as 
a source material since it simultaneously (< 10 psec) emits a v-ray along 
" 1th the pOSltron - This 1*2B HeV T-ray acts to start the timing apparatus 



which is stopped by the arrival of the - 511 KeV y-ray signifying the occur- 
rence of the annihilation event. The positron lifetime depends on the con- 
centration and distribution of the electrons in the region of the material 
through which the thermalized positron drifts before Its annihilation. This 
"mean- lifetime will of course be affected by any process, such as the 
trapping of positrons at defects, which alters the local concentration of 
electrons seen by the positron (Coleman 1977). These mean lifetimes, which 
are on the order of hundreds of picoseconds, can be measured to an accuracy of 
a few picoseconds with modern equipment. 

As mentioned previously, when the annihilation event does occur the 
electron-positron pair will usually have a small net momentum in the labora- 
tory frame. This net momentum affects both the energy of the resulting Y -ray 
as well as its direction of emission as shown In Figure 1. The relationships 
between the momentum of the annihilating y-ray, the measured Doppler shift in 
energy, and the measured angle 5* are described in Appendix B. 

\ Center of Mass 
\ Velocity V 




gure 1. A Schematic Diagram Showing the Effects of the Electron-Positron 
Pair Momentum on the Energy and Oi recti on of the Two Resultinq 
Annihilation y-Rays in the Laboratory Reference Frame 
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Angular correlation measurements, which measure the deviation from co- 
linearity of the two photons emitted from the annihilation process, were first 
reported by DeBenedetti and his coworkers in I960. The so called long-slit 
geometry (DeBenedetti 1950) is commonly chosen since it provides both good 
angular resolution and adequate counting rates. However, other geometries can 
be employed for specific applications. The angular distribution of a sample 
is determined by measuring the coincidence counting rate as a function of 
angular deviation from colinearity. In well annealed metals this data gener- 
ally takes the form of a central Inverted parabola with Gaussian wings. The 
Instrument resolution is good, typically ~ 0.5 mrad, but the counting rate is 
comparatively slow and the required source strengths are relatively high. 

Alternatively, Ooppler broadening measurements could be used to deter- 
mine the net momentum distribution of the annihilating electron -positron pairs 
by measuring the energies of the resulting Y-rays. The center of mass momen- 
tum of the annihilating pair along the line of emission of the Y-ray gives 
rise to Doppler shifts of the Y-ray energies. The distribution of these 
momentum components then results in a Doppler broadening to give the 511 KeV 
line a width of approximately 511 (V F /2c) KeV where V F (- 108 cm / s ) i s the Fermi 
velocity of the conduction electrons and c 1s the speed of light (Hughes 
1982). Details of the experimental procedure for collecting Doppler-broadened 
spectra are given in Chapter II. 

Although lifetime measurements, angular correlation, and Doppler- 
broadening measurements are by far the most common PAS techniques there are 
two others which should be mentioned. One of the drawbacks of using conven- 
tional sources such as 22 Na or 68 Ge is that only certain depths of the mater- 
ial under investigation are sampled extensively by these positrons. Hughes et 
_.„ r ^ 4lvJpa2X^ ... . „ 
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« - *J-M 0 C 2 (6e) = t 255.5 («6)ev (K2) 
where 6E is the Doppler Energy Resolution in eV 
and 6B is the Angular Correlation Resolution in mrad. 

The best Doppler systems have a resolution of about 1.2 KeV which 
corresponds to an angular correlation resolution of 5 mrad. Typical angular 
correlation systems have a resolution approximately an order of magnitude 
better than this. For Doppler broadening techniques to effectively compete 
with angular correlation techniques in high precision applications it is 
necessary to improve the effective resolution of Doppler PAS. Schultz et aK 
OKI) point out that theoretically a 10-fold increase In effective resolving 
power of Doppler measurements can be achieved by deconvolution. 

Deconvolution is the mathematical technique used to separate the 
instrumental broadening effects from the true Doppler broadened positron 
annihilation Y -ray spectrum. Instrumental broadening is a complicated problem 
caused by many different factors. For example, physical properties of the 
detector crystal such as frontal dead layers, geometry of the active region 
and regions of poor charge collection, can all lead to line broadening. The 
Primary contribution to instrumental broadening in a solid state semiconductor 
detector is due to statistical fluctuations in the division of absorbed energy 
between ionization and heating of the crystal lattice, i.e. the Fanno effect. 
(McNelles and Campbell 1975). 

The shape of the instrument resolution function, as well as the mag- 
nitude of the line broadening, is most often determined from the 514 KeV Y - ra y 
emitted from the isotope 8S Sr . The Mffiai ^ „ M f|< ^ ^ ^ ^ 

Ooppler-broadening of the 511 KeV Y - ra y line in the absence of defects in a 
metallic crysta! is_of the same order <~ KeV) as the detector resolution, and 
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changes Induced by defects are significantly smaller. Fortunately, the defect 
induced changes in the deconvoluted Y-ray energy distribution are much more 
dramatic than 1n the instrumental ly smeared raw data (McNelles and Campbell 
197b). In addition, removing instrumental effects offers another advantage. 
Such a procedure allows direct comparison of experimental results obtained in 
different laboratories. If results are reported without deconvolutlon they 
must be nondimensionallzed in some way or else the Instrument resolution func- 
tion must also be reported. In either case the interpretation of results is 
complicated in the absence of a deconvolutlon procedure. 

In light of the above discussion it was decided that the results pre- 
sented 1n this work would be for deconvoluted Doppler spectra. Unfortunately, 
the mathematics involved are somewhat complicated. Because of this, most of 
the existing deconvolutlon schemes are iterative in nature and hence computa- 
tionally inefficient. One of the major goals of this work was to develop a 
deconvolutlon scheme for Doppler broadened positron annihilation measurements 
which is model Independent, does not require iteration, and which yields 
reproducible results. 

In addition to the need for deconvolutlon, OBPAS also currently suf- 
fers from an Inability to distinguish various defect types. For example, 
vacancies and dislocations tend to alter the annihilation lineshape in quali- 
tatively the same manner. That is, in both cases the spectrum becomes more 
narrow as the defect concentration increases. Most of the lineshape param- 
eters currently used to characterize Doppler broadened spectra are unable to 
distinguish the effects of different defect types because they describe only 
gross changes in lineshape. For example, Mackenzie's S-parameter, which is 
the most commonly used lineshape parameter, essentially measures the width of 
the distribution (Mackenzie et al. 1970). Hence, for an V ni««,n c^i.™ ,..-u 
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an analysis yields a single value. This is useful if there is only one degree 
of freedom in the series of experiments. If two variables are allowed to 
change, however, it is impossible to separate their effects* Clearly, 
additional parameters are required. The other parameters selected must be 
uncor related with the S-parameter. That is, if the second parameter chosen is 
a function of the first parameter only, then no new information is gained. 
Campbell (1977) and Schultz et aU (1981) have discussed the use of other 
possible parameters to describe Uneshapes. Although these parameters may be 
more sensitive to variations in specific defect concentrations they are 
correlated with S, and therefore of little help in differentiating defect 
types • 

In light of the above discussion it can be seen that it would be de- 
sirable to find a family of lineshape parameters which is capable of di sting- 
ui shiny various types of metallurgical defects. For example, one of the 
groups of parameters to be examined in this work Is the central moments of the 
Y-ray energy distribution. It 1s well known that the values of these moments 
give Information regarding the shape of the distribution. In fact, it can be 
shown that if all of the moments are known one can construct the distribution. 
The even moments give information about the width of the distribution while 
the odd moments vanish for symmetrical distributions such as the ones to be 
examined in this work (the symmetry question Is addressed in Chapter III). 
The second central moment is called the variance and its square root is the 
familiar standard deviation. The fourth central moment is called the Kurtosis 
or flatness. The advantage of moment analysis is that the information ob- 
tained is not redundant^ For example, it is very possible to have two dis- 
tributions with equal standard deviations (second moments) and yet widely 
di fferent fourth moment s. Such a family of lineshape parameters suggests the 
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possibility of separating the effects of multiple metallurgical defect types 
using OBPAS. 

The above brief review of positron annihilation spectroscopy as an NOT 
tool for the study of metallurgical defects indicates that the technique with 
the greatest potential for practical use outside the laboratory 1s uoppler- 
broadenlny. This technique offers the advantages of low source strength re- 
quirements, the ability to use 6 «Ge as a source material In order to reduce 
surface effects, and few geometrical experimental constraints. Doppler broad- 
enlng currently suffers from relatively poor energy resolution and an Inabil- 
ity to distinguish different metallurgical defect types. Based on these ob- 
servations, the purpose of this work is twofold: 1) to develop an efficient, 
model independent, noniteratlve deconvolution method, and 2) to develop a 
family of Hneshape parameters capable of distinguishing various defect 
types. 



CHAPTER II 
EXPERIMENTAL PROCEDURES 



This chapter Is divided into two sections. The first section des- 
cribes the Doppler broadened PAS instrumentation and the methodology employed 
to collect Doppler spectra. The second section outlines the thermomechanical 
treatments used to create the desired defect concentrations in the metallic 
samples which were studied in this work. 

The Popper PAS System 

The heart of the Doppler broadened PAS instrumentation system is the 
single crystal high purity Germanium (HPGe) detector. This detector is mount- 
ed at the end of a cold finger in a cryostat and, to reduce thermal noise, is 
operated at liquid nitrogen temperatures. The cylindrical detector crystal 
used in this work is 46.7 mm in diameter and 46.2 mm 1n length. When the 
annihilation y-rays impinge upon the crystal their energy is absorbed by two 
mechanisms. The first is the creation of electron-hole pairs within the 
crystal. These electron-hole pairs are subsequently swept apart by the 
applied detector bias (4000 V) and flow through an external circuit creating a 
current pulse whose integral with respect to time has an amplitude propor- 
tional to the incident Y-ray energy. The second absorption mechanism is heat- 
ing of the crystal. The ratio of these two absorption mechanisms is known as 
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the Fanno factor and the statistical fluctuations in the division of absorbed 
eneryy between ionization and heating is the primary contribution to instru- 
mental line broadening in solid state semiconductor detectors. For further 
details of the absorption process see Moore (1983). 

The resulting current pulse from the detector is rather weak and un- 
suitable for cable transmission. Because of this a preamplifier (Ortec Model 
120-5) is mounted directly adjacent to the detector output. The detector and 
preamp are U.C. coupled, as is the entire system, and powered from the same 
A.C. source to prevent ground loops. The preamp is a charge sensitive device 
capable of accepting signals corresponding to the y-ray energy range of 0-10 
MeV. Its output is a voltage pulse with an amplitude proportional to the 
integrated charge supplied by the detector and hence proportional to the Inci- 
dent Y-ray energy. The magnitude of the voltage pulse amplitude Is related to 
the incident y-ray energy by the approximate conversion factor of 280 mV per 
MeV. Therefore for a typical 511 KeV y-ray the output of the preamp would be 
approximately 143 mV. In addition to converting the detector charge pulse to 
a voltage pulse with an intensity suitable for cable transmission the preamp 
also shapes the resulting voltage pulse so that 1t is acceptable to the linear 
amplifier. 

The linear amplifier used in this work is an Ortec Model 571. It 
accepts the shaped pulses from the preamp and amplifies them to suitable volt- 
ages for transmission to the multichannel analyzer. The output pulses are 
active filter shaped (skewed Gaussian approximations) and have amplitudes in 
the 0-10 V range. Unipolar outputs are used to optimize the signal-to-noise 
ratio (resolution). Since the main amplifier is operated within its linear 
range the amplitude of the output voltage pulse is proportional to the inci- 
dent y-ray en ergy as measured b y the HPGe detector. 
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The output from the linear amplifier is routed to an Ortec Model 408A 
biased amplifier. This device accepts the output pulses from the main ampli- 
fier and allows one to selectively view a:.y region of the Initial 0-10 MeV en- 
ergy range. In this work the biased amplifier was used to select the energy 
window of 483 to 548 KeV. The output pulses from the biased amplifier serve 
as the Input to the analog-to-d1gital (ADC) converter located within the Nor- 
land IT-5300 multi-channel analyzer (MCA). In this device the amplitude of 
the voltage pulse 1s converted to a digital representation of the energy of 
the detected Y-ray and subsequently dynamically displayed on the MCA's CRT in 
the form of a histogram of energy versus count intensity. The MCA can be pro- 
grammed to accept input signals until a predetermined number of counts 1n the 
peak channel is reached. At this time the experiment is terminated and a com- 
plete Ooppler broadened spectrum Is displayed on the CRT. The energy distri- 
bution is then ready for transmission to the computer for subsequent data man- 
ipulation. 

The two most Important parameters which characterize the quality of a 
Ooppler PAS system are the width, in units of energy, of the Instrument reso- 
lution function and the energy resolution, in units of eV per channel, of the 
MCA. The width of the instrument resolution function (IRF) is, in essence, a 
measure of the quality of the detector crystal, although 1t 1s also affected 
by the associated amplification electronics as well as the ADC. The IRF is 
most often determined by examining the response of the system to the monoener- 
getic (E = 514 KeV) photons emitted from the isotope 85 Sr. If the instrumen- 
tation was perfect the energy distribution, as recorded by the MCA, would con- 
sist of a single bar in the histogram in the channel corresponding to 514 KeV. 
However, in any real instrumentation this line will be broadened for reasons 
P. rey ig^- glassed y-- Th - e , am ount . of in^rumental broajteMnu is mgsl of ten,™- 
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ported in terms of the energy width of the 85 Sr line as measured across the 
two channels which contain half the number of counts in the peak channel. 
This energy measure Is commonly named the full width half maximum (FWHM). 
Figure 2 shows typical IRF for the Doppler PAS system used in this work. As 
can be seen in this figure the FWHM of the IRF Is about 1.28 KeV. 

The second parameter which describes the quality of the Doppler PAS 
instrumentation 1s the energy resolution of the MCA which is determined by the 
following procedure. First, a Doppler broadened spectrum is collected and the 
channel corresponding to its centroid is determined. The centroid channel is 
assigned the energy value 511 KeV. Next, an energy spectrum 1s obtained for 
the Isotope 85 Sr. The channel corresponding to the centroid of this distribu- 
tion is calculated and assigned the energy value 514 KeV. The MCA energy 
resolution is then easily calculated by dividing the energy difference (3 KeV) 
by the number of channels between the centroids. In this work the 511 KeV 
channel was number 439 while channel number 486 corresponded to 514 KeV. 
Hence, the MCA energy resolution Is 3000 eV/(486 -439) or - 64 eV per channel. 
Clearly this value can be adjusted using the biased amplifier. The reason for 
selecting the 64 eV/channel is to obtain the best possible resolution while 
maintaining the ability to display the entire region of interest (483 to 548 
KeV) on the 1024 channel MCA. The reasons for requiring the above energy 
window will be discussed in Chapter III. 

While the MCA energy resolution is a parameter which, to a certain ex- 
tent, can be varied, the IRF is essentially predetermined by the instrumenta- 
tion. The one factor contributing to the width of the IRF which can be con- 
trolled is the data acquisition, or count rate. It is well known that in- 
creasing the count rate decreases the instrument resolution. Therefore, one 
m " St * e l e £ t ., ? w hich . results i n an acc eptable line broa den i n 9. w hile 
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keeping the time required to collect a spectrum reasonably short. Moore 
(1983) demonstrated that with the Doppler instrumentation used in this work 
very reproducible results can be obtained with count rates of < 2 kHz. The 
count rate Is easily adjusted by changing the sample to detector distance. 
All of the results reported in this work were obtained using count rates of 
< 2 kHz with an associated experiment time of ~ 30 minutes resulting 1n a peak 
height of 40K counts and - 3 x 10* total integrated counts. 

The positron source isotope in all experiments was 68 6e which was 
chosen to avoid surface effects. The mean depth of penetration in aluminum 
for such positrons is ~ 0.50 mm. The activity of the source was initially 50 
uC. The half-life of 68 Ge is _ 275 days and hence ^ activUy Qf ^ source 

did not change significantly for any single set of experiments since a typical 
set of measurements requires less than a week to accumulate. It should be 
pointed out that the activity of the source does not measurably affect the 
observed annihilation spectra, provided the count rate remains constant. 

Although the count rate can be easily adjusted by changing the source 
to sample distance such a change introduces an additional error source. This 
error source Is related to the source/sample geometry. It has been demon- 
strated that the observed spectra may depend on the angle with which the gamma 
rays impinge upon the detector crystal (Smedskjaer 1983). It has also been 
shown that In high precision experiments among samples of the same material, 
one should use samples of nearly identical thickness to reduce the error asso- 
ciated with small angle scattering (Smedskjaer 1983). 

The 68 Ge source was deposited on a thin (~ 0.05 mm) nickel foil for 
handling convenience. This results In what is known as a source contribution 
to the annihilation lineshape. That is, in each experiment some of the posl- 
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events will contribute to the overall annihilation spectrum. However, the 
ratio of source annihilation events to sample annihilation events Is small (< 
1%), and the source contribution is neglected in the remainder of this work. 

In all measurements the source and sample were arranged in what is 
commonly called the sandwich configuration. In this configuration the posi- 
tron source is placed between two identically treated metallic samples. The 
samples are thick enough, typically a few millimeters, to guarantee that 
virtually all of the positrons annihilate within the material being examined. 

Specimen Preparation 

The material examined in this work is high purity, large grain size 
aluminum purchased from the Materials Research Corporation. The aluminum was 
zone refined to a nominal purity of 99.9995% before being vacuum melted and 
cold fabricated to size. The chemical analysis, in PPM of metallic inter- 
stitials, is listed in Appendix C. The aluminum was supplied as a 13 ran 
diameter rod from which ~ 1.5 mm thick disks were machined. Each of these 
disks had an ~ 1.5 rnn hole drilled through it near its edge. The purpose of 
this hole was to allow the disks to be suspended from a titanium wire while 
within the vertical tube furnace. After the disks were machined to size, they 
were given a full anneal heat treatment consisting of 2 hours at 400°C fol- 
lowed by furnace cooling overnight. The as annealed aluminum disks were the 
starting material for all subsequent thermal treatments. 

The thermal treatments used in this work are very similar to those 
outlined by Cotterill et al_^ (1972). In this work the authors describe a two 
step time/temperature cycle which allows one to produce, in a controlled 
fashion, first samples > witj^ vacancies as the domi nate defect typ_e*_ and_tbea t 



samples with dislocation loops as the major defect type. The actual thermal 
treatment employed in this work is as follows: 

1) Heat the samples to the desired quenching temperature and hold 
them at temperature for one hour to establish thermal equilibrium. 

2) Quench the samples from the furnace into - 293°K water. 

3) Ory the samples within 5 to li seconds. 

4) Quench from - 293°K to a liquid nitrogen bath at 77°K. This 
procedure has been shown to produce primarily vacancies, the concentration of 
which 1s a direct function of the quenching temperature. 

5) After DBPAS characterization, the aluminum disks are then annealed 
at 3b3°K for thirty minutes, which results in the production of prismatic 
dislocation loops by vacancy coalescence. 

6) The samples are cooled to 77°K and the PAS measurements corres- 
ponding to the dislocation loops are made. 

7) Finally the disks are annealed at S23°K for two hours to remove 
the dislocations and return the material to the annealed state. 

Several of the steps in this procedure require further explanation. 
The quenching procedure Is taken from Cotterill et aU (1972) and their 
results together with the following discussion forms the basis for the 
assumption used in this work that the PAS measurements made after step 4 are 
for a fairly uniform distribution of vacancies of concentration approximately 
equal to the equilibrium concentration at the quenching temperature. 

Quenching experiments are usually performed on thin metallic foils so 
that the quench rate remains constant across the sample. However, as dis- 
cussed previously, constraints imposed by the positron annihilation experi- 
ments put a lower limit on the sample thickness of ~ 1.5 mm. For a complete 
review of the complications which can evolve during uuenchino t n 
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sample thickness see, for example, Cotterill et art. (1965). The main problems 
to be addressed here are the formation of vacancy clusters and the loss of 
vacancies at the sample surfaces. The second problem was Investigated in a 
paper by Wilkes and Harrison (1972) in which they observed the quenching of 3 
mm thick disks of aluminum and reported a very low vacancy loss. Since the 
samples used in this work are only ~ 1.5 mm thick and since the quench pro- 
cedure is similar to that employed by Wilkes and Harrison, 1t will be assumed 
that the vacancy loss during the quench can be ignored. As for the formation 
of vacancy clusters, the positron lifetime measurements of Cotterill et al_. 
(1972) demonstrate that for this quenching procedure the probability of form- 
iny bivacancies, or larger clusters, is small. 

It is very important to dry the samples thoroughly between the water 
quench and the liquid nitrogen quench. This 1s necessary to avoid the forma- 
tion of a thin layer of ice on the specimen surface which would contribute an 
undesirable component to the measured PAS spectrum. As pointed out by 
Cotterill et aU (1972) the annealing of the vacancies during this short dry- 
ing procedure may be neglected since the mean lifetime of vacancies (at room 
temperature) in a specimen quenched from 873°K 1s ~ 100s and increases for 
lower quenching temperatures. It 1s, of course, desirable to keep the drying 
time as short as possible. 

Step 5 of the procedure, the formation of prismatic dislocation loops 
from the collapse of disk shaped vacancy clusters during a low temperature 
anneal of aluminum, was reported by Hirsch et a]_. in 1958. Cotterill et al. 
(1972) demonstrate that the appropriate thermal treatment consisted of a 30 
minute anneal at 353°K. The theoretical relationship between the vacancy 
concentration after the quench, C v , and the concentration of dislocation 
cor es, Cq- (gi ven as an a tomic fraction), is presented in Appendix G. 
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The apparatus used to perform the thermal treatments is shown In 
Figure 3. A quartz tube passes through a vertical tube furnace and is 
attached to the room temperature water reservoir by a "quick disconnect" 
coupling. Just above the water level in the reservoir Is the argon supply 
port. All thermal treatments were performed in an argon atmosphere to avoid 
oxidation. The other end of the tube is connected to a flow meter which 1s 
used to monitor the argon flow rate. The other side of the flow meter is 
vented to atmosphere. The temperature is measured by a chromel-alumel thermo- 
couple positioned Inside the tube furnace and adjacent to the aluminum disks. 
These disks are supported by a titanium wire, the other end of which is con- 
nected to an iron slug. The slug, and hence the specimens, are held in posi- 
tion in the furnace by a magnet located outside of the tube. The actual 
quench 1s performed simply by pulling the magnet away from the tube and allow- 
ing the entire specimen-wi re-slug arrangement to fall by gravity into the 
quench bath. It Is, of course, important for the system to be designed so 
that the specimens fall freely through the tube and quickly reach the quench 
bath. If the specimen-w1 re-slug arrangement does not fall freely the tempera- 
ture of the samples when they enter the water will be significantly lower than 
the temperature measured by the thermocouple. 
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CHAPTER III 
DECONVOLUTION OF DOPPLER PAS MEASUREMENTS 



As discussed 1n the Introduction, compared to angular correlation (AC) 
measurements, the major disadvantage of the Doppler PAS technique is its in- 
herently poor energy resolution. In order for Doppler methods to be competi- 
tive with AC methods in high resolution applications the instrumental line 
broadening must be removed from the experimentally obtained Doppler spectra. 
Mathematically, the technique which accomplishes this separation is deconvolu- 
tion. In this chapter the deconvolution of Doppler PAS measurements will be 
discussed. 

The chapter Is divided into four sections. In the first section the 
mathematics of deconvolution are presented. This section justifies the use of 
deconvolution techniques to separate instrumental effects from the desired 
true Ooppler-broadened energy spectrum. In addition, the inherent advantages 
and mathematical difficulties associated with the use of Fourier analysis to 
solve the problem are addressed. 

In the second section of this chapter the methods introduced by other 
investigators to solve the incorrectly posed deconvolution problem, along with 
their inherent advantages and disadvantages, are discussed. The third section 
presents the FFT/power spectrum deconvolution method (Schaf fer et aj_. 1984); 
and describes its use in detail. In the final section of this chapter the 
FFT/power spectrum method is demonstrated using a neutron irradiated high 
purity single crystal sapphire sample. 
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The Mathematics of DeconvoTution 

The first objective of this section Is to demonstrate the validity of 
using the convolution equation to describe the relationship between h(x), the 
experimentally observed energy distribution, f(x), the true Doppler broadened 
energy distribution, and r(x), the instrument resolution function { IRF) . In 
Figure 4(a) the instrument resolution function r(x) is shown. Recall that this 
function represents the instrumental broadening of what should be a single 
energy line at the centroid of r(x). The function r(x) is assumed to be norm- 
alized to unit area. In order to understand the relationship between these 
three functions, f(x) should be viewed as a histogram comprised of segments of 
width Ax. If it can be assumed that the instrumental response is linear (the 
detector and all associated electronics are operated within their linear range) 
then each segment of the f(x) histogram will be broadened by the unit area IRF 
r(x). Hence, the value of the observed function h(x) at x 1 will have contribu- 
tions from many of the f(x) histogram segments in the neighborhood of x 1 . Fig- 
ure 4(b) is a schematic representation of the contribution to h(x') from the 
histogram section of f(x) at a distance xi from the origin. This contribution 
Ah^x 1 ) will be the product of the area f (xj )Ax and rCx'-x^, that is: 



The total value of h(x') will be the sum of the contributions from all pos- 
sible histogram segments in f(x). Thus, 



or as the width of the segments approaches zero the summation in Eq. 3.2 is 
replaced by an integral to yield, 



Ahi(x') » f(x 1 )r(x' - x,)Ax 



(3.1) 



b(x') = I Ah^x 1 ) = I ffx^rU 1 - Xi )Ax 




(3-2) 



CO 



H(x') = / f(x)r(x' - x)dx 



(3.3) 
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Clearly the same procedure holds for any value of x In h(x) and the final form 
of the relationship, with a change of variables, is: 

cm 

h(x) = /. - y)dy (3>4) 

Eq. 3.4 is the convolution equation which serves as the starting point for the 
remainder of the mathematics used in this work. 

The advantage of using Fourier analysis when dealing with convolution 
problems can be seen clearly when Eq. 3.4 is rewritten as: 

H(w) - F(w).R(w) (3>5) 
where the functions In Eq. 3.5 are the Fourier transforms of the respective 
functions in Eq. 3.4 according to the general relationship: 

■M £ e-i" 9 (x>dx (3 . 6) 

Note that in Fourier space convolution is simply a complex multiplication. 
The equivalence of Eqs. 3.4 and 3.5 is demonstrated in Appendix 0. 
Rearranging Eq. 3.5 yields 

F(w) = H(w)/R(w) (3#7) 
from which one can find f(x) using the Inverse Fourier transform 



Eq. 3.8 is an expression for the true Ooppler broadened energy spectrum In 
terms of the Fourier transforms of the experimentally observed function h(x) 
and the instrument resolution function r(x). However, the uniqueness of the 
solution given in Eq. 3.8 requires that h(x) be known exactly. This condition 
1s never satisfied in a real experiment and this fact significantly compli- 
cates the mathematics of the problem. The observed function h*(x) will be 
related to the "true" function h(x), as defined in £a. 3.4. fnr a „ „„, 
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function f(x) by: 

max |h*(x)-h(x)|<e (3.9) 
-«<y<+» 

where e is a measure of the noise associated with the stochastic nature of the 
experiment. 

It can then be shown (see Appendix E) that there are two solutions 
f](x) and f2(x) ■ f](x) + C s1n(2«nx/(b-a)) with corresponding function h^ and 
h2 (as defined by Eq. 3.4) both of which satisfy the inequality (3.9). Hence, 
both f] and f£ must be considered as solutions to Eq. 3.4 with a left hand 
side h*(x) taken from an experiment to within an error e (Turchin et al . 
1971). In fact, there is an infinite set of such solutions a subset of which 
contain indefinitely large rapidly oscillating components. 

In order to understand the mathematical source of these difficulties 
Eq. 3.8 was investigated further. If in Eq. 3.8 H(w) is replaced by H*(w) 

with <• - • 

H*(")« — J e lwx h*(x)dx (3.10) 

and h*(x) = h(x) + 6(x) (3.11) 

where 6 is a stationary random process with zero mean, then it has been shown 
(Kautian 1958) that the dispersion of f(x) is 

00 

D(f) = <f 2 > - <f> 2 = — J ^"K dw (3.12) 

m — iR(w>r 

where G(w) is the power spectrum of 6(x) 

and |R(w)| + 0 for |w| ♦ - 
Therefore the dispersion of the solution f(x) can remain finite only if the 
power spectrum of the noise approaches zero rapidly as w approaches infinity. 
In practice this will not be the case since the noise will inevitably contain a 
"white noise" component such that G(w) remains finite as |w| + * and the dis- 
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persion of the solution f(x) therefore becomes infinite (Turchin et aj_. 1971). 

Alternatively one can understand the source of the mathematical diffi- 
culties as follows. The observed function h*(x) contains high frequency noise 
components which are not physically realizable in the true function h(x). In 
Eq. 3.8 these high frequency components are divided by small eigenvalues, R(w), 
which results in large oscillations in the solution f(x) (Turchin et aU 1971). 

Recognizing these difficulties, previous investigators have proposed 
several deconvolution techniques which appear to successfully solve the prob- 
lem. Since the accuracy of the deconvoluted spectrum can never really be de- 
termined, except in certain carefully constructed numerical experiments, the 
criticisms of a specific method tend to relate to the amount of computer time 
required, the complexity of the programming effort, the degree to which the 
analysis may be automated, or to the need for some a priori information con- 
cerning the function f(x). As a preface to the introduction of the FFT/power 
spectrum deconvolution technique, the following section 1s a review of the 
prior techniques that have been developed to deconvolute Ooppler PAS data. 



Previous Deconvolution Techniques 

One of the early Ooppler PAS deconvolution schemes Introduced by Hotz 
et aK (1968) is basically a matrix approach. The convolution in its discrete 
form 1s 

a 

b(x) - I r(y - x)f(y)Ay { 3j 3) 
y=-a 7 

where Ay is the incremental width (resolution) 

and the other functions have their usual meanings 

This can be rewritten in matrix form as 
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b i * A j1 c j (3.14) 
where b is the discrete experimental vector (spectrum) 

c is the unknown discrete vector before instrument broadening 
and A is a rectangular matrix 1n which each row Is the detector re- 
sponse function appropriate to the corresponding component of c. 
in the method developed by ttotz et aK (1968) the matrix A is determined using 
«5 Sr and the vector b is the output of the Doppler system MCA for a given 
experiment. In this technique an Initial estimate of the unknown vector c is 
made, say cO), and a corresponding vector b<D is calculated from Eq. 3.14. 
The quality of c(D is determined by finding X 2 for the fit between b and b<D 
as defined by 

X 2 = X(bO) - bO/*>i (3.1b) 

Next cd) is modified in a manner described in the original work to create a 
second estimate, cU), of the desired function c. Then bU) 1s found from Eq. 
3.14 and a new X 2 value determined. The process Is repeated until X 2 is mini- 
mized. This process is extremely involved and inefficient. In addition the 
Initial estimate cO) as well as the modification criteria used to generate 
subsequent c«) estimates can influence the final result. 

In more recent years the matrix approach has been examined in greater 
detail (Aulick and Gallie 1981). Although this work 1s net related to Doppler 
PAS measurements, the mathematical problem is the same. It 1s demonstrated 
that 3.14 should be rewritten as: 

b=Ac (3.16) 

where b = b + e 

and c is the measurement error 
The least squares soluti on of Eg. 3. 16 is then 



on 
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Atb ■ c (3.17) 
where A* is the pseudo-inverse (Penrose 1956) of A. 
However, due to the fact that b contains noise Eq. 3.16 is incorrectly posed 
and A will be ill-conditioned. Hence, the elements of At are large and the 
error in b will be magnified, obscuring the true solution which Is 

c * Atb (3.18) 
Aulick and Gallie (1981) point out that to avoid this problem one needs to 
calculate solutions to Eq. 3.16 in a way which will damp out the error in b. 
One way of accomplishing this goal Is to supply some additional a priori in- 
formation about the shape of the desired vector c. The additional Informal 1 
recasts the problem as a constrained optimization problem. This technique is 
called regularization. 

An example of the type of a priori information commonly supplied 
appears in the deconvolutlon technique used by Reddy and Carrlgan (1970). In 
their approach the resolution function is assumed to be the sum of two 
Gaussians with the same centroid but differing In their amplitudes and widths. 
The shape of f(x) (returning to the more common notation) is assumed to be a 
truncated parabola based on the free electron model (Reddy and Carrlgan 1970). 

Introducing these constraints on the shapes of r(x) and f(x) into the 
least square solution problem yields a physically plausible result. That is, 
the function f(x) which is obtained by this method does not contain indefi- 
nitely large and rapidly oscillating components typical of the results obtain- 
ed by evaluating Eq. 3.8 directly. The difficulties associated with this 
method are related to the validity of the a priori assumptions. For example, 
the resolution function r(x) is constrained to be symmetric. In addition as- 
suming the shape of the true Ooppler energy distribution, f(x), to be a simple 
* ™"" ted Parabola, j s. pjobably an ov ers.1«p llf-tcation»- Qther^AftvjesMg 
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(MacKenzle et aj_. 1976) suggest the modeling of f(x) using the superposition 
of Gaussian and inverted parabolic functions as representing contributions 
from core and valence electrons respectively. While it is not the purpose of 
this work to discuss the validity of the various electron distribution models, 
the very fact that there exists more than one analytic prediction of the shape 
of f(x) Is strong reason to develop a deconvolution technique which places 
minimal constraints on the shape of f(x). 

One of the most widely used deconvolution techniques in Doppler PAS 
studies is an iteratively based method Introduced by Jackman et a]_. (1974). 
In this method an Initial estimate of the desired result f(x), say f^x), Is 
selected. The estimate fj(x) is then convoluted with r(x) to obtain h^x), 
the first estimate of the experimentally obtained energy distribution h(x). 
The function hj(x) is then compared point-by-point with h(x) and if satisfac- 
tory agreement is obtained fj(x) is assumed to be a reasonable approximation 
to the true solution. If unsatisfactory agreement Is obtained then fi(x) is 
modified to produce f 2 (x) and the process is repeated. 

The initial estimate of f(x), the modification procedure for obtaining 
f 2 (x) from f^x), and the criteria for stopping the iteration procedure are 
all current topics of investigation 1n the literature. In the original method 
the experimentally obtained function h(x) served as the initial estimate 
f](x). The modification algorithm was 

*2(*1> " fi(xi) x ( T h(xj)/ T h^x,)) (3.19) 

j=i-2 J j=1-2 J 

Subsequent work demonstrated that the convergence rate of the Iteration pro- 
cess could be improved If the summations in Eq. 3.19 were taken over nine 
channels instead of the original five and if the summation ratio was squared 

** f0r *,™}* :£] 1 cat : j™ frr- lliH) (^Mgfae^QgJgCCUaZ^ THp iter^ijn^ro^^ 
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cedure is stopped when the comparison of h(x) and h N (x) yields a xZ value less 
than some predetermined minimum (Jackman et ^1_. 1974) or when the sum of the 
squares of the deviations falls below some predetermined level (Dannefaer and 
Kerr 197S). 

The main advantage of these iteratively based techniques is that a so- 
lution is obtained without any restriction on the shape of f(x). While this 
1s indeed an attractive feature of the method it also suffers from some inher- 
ent disadvantages. These disadvantages Include the amount of computer time 
required, the possible dependence of the result on the initial estimate of 
f(x), and the possibility of introducing large oscillations Into the solution 
if the stopping criterion is selected with too small an error bound. This 
last difficulty Is of course related to the mathematical nature of the incor- 
rectly posed problem. 

The iterative deconvolution technique essentially Ignores the statis- 
tical nature of the problem. An alternate approach is to treat Eq. 3.4 as a 
problem of mathematical statistics and introduce a regularizer (a priori In- 
formation) in the form of some probability distribution for f(x) (Turchin et 
aU 1971). This 1s the method of statistical regularizatlon which offers the 
advantage of error analysis. The result of this process Is the most probable 
f(x) for a given h(x) measured with mean square errors S(x) which are indepen- 
dent of x and normally distributed with zero mean (Goldanskii et al_. 1978). 
The regularizer 1s a form of a priori information concerning the shape of f(x) 
and hence somewhat restrictive. However, such an approach appears to yield 
reasonable results if a solution can be obtained using minimal a priori infor- 
mation concerning the smoothness of f(x). 

In 1981 the Fourier cut-integral method was introduced by Szotek and 
C J}$?. e lewsklt r In this •-techntqu^ l Eq ;<i --3-*AdA = cewr1 ttenr-as 
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The difference in the equations is the introduction of finite limits of inte- 
gration. If p 0 Is then chosen small enough the solutions obtained are physi- 
cally plausible. The mathematical explanation for this result Is that the in- 
sertion of finite limits of Integration forces the dispersion of f(x) to 
remain finite. (Refer to Eq. 3.12 and the associated discussion). The prob- 
lem associated with this method is that the optimum value of p 0 is not readily 
predictable. Approximations for f(x) are obtained using different values for 
p 0 . and one must then select the "best" approximation to the true f(x). This 
is typically done by convoluting the function with r(x) and then comparing the 
result with h(x). As p 0 becomes too large the solution begins to exhibit 
large oscillations as would be expected from the mathematical theory. 

In 1982 Verkerk published an extensive mathematical treatment of the 
deconvolution problem. In this work he discusses the relative strengths and 
weaknesses of many of the currently employed deconvolution techniques. Ver- 
kerk then proposes a method which in many ways is very similar to the Fourier 
cut-Integral method. The crux^of this method is that the high frequency 
Fourier coefficients of h(x) do not contain significant information and conse- 
quently only a limited number (say K) of Fourier coefficients H(w) contain the 
essential information from measurement of h(x). It should be pointed out that 
if H(w) contains N Fourier coefficients then deleting the N-K high frequency 
coefficients in the reconstruction of f(x) is equivalent to Inserting finite 
limits of integration in Eq. 3.4 as in the Fourier cut-integral method. 

The deletion of the high frequency coefficients is equivalent to a 
rectangular low-pass filter and it is well known that the application of such 
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a filter to F(w) will cause nonphyslcal oscillations in the approximation to 
f(x). 

In summary, there are two different types of errors associated with 
the estimate of f(x): 1) statistical errors, which are amplified by the de- 
convolution operation and 2) systematic errors due to filtering, which are 
largest in the peak region and may cause nonphyslcal oscillations in the tails 
of the deconvoluted signal. Error type 1 increases as the number of Fourier 
coefficients used in reconstruction increases, while error type 2 decreases 
(Verkerk, 1982). Thus, by using such an approach the problem is reduced to 
selecting the optimum number of coefficients. One disadvantage of Verkerk' s 
method, as with the Fourier cut-integral method, 1s that the procedure must be 
repeated several times for different numbers of coefficients until the optimum 
number Is determined. It would be desirable to develop a simple one-step 
method to determine the appropriate number of coefficients to use in recon- 
structing f(x) from F(w). and indeed this is one of the primary advantages of 
the FFT/power spectrum deconvolution method which 1s described in the follow- 
ing section. 

The FFT /Power Spectrum Deconvolution Method 

The deconvolution method developed 1n this work is a fast Fourier 
transform (FFT) based method which is simple to execute, fast, easily auto- 
mated, and assumes little or no a priori information concerning the shape of 
f(x). The distinctive feature of this method is that the number of Fourier 
coefficients retained in the reconstruction of f(x) is determined by examining 
the power spectrum of the experimental data. Using the observation that the 
power spectrum falls steeply to a minimum before rising due to the white noise 
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coefficients at this first minimum. The logic behind this approach is rather 
simple. Retaining fewer coefficients than this implies a smoothing of the 
deconvolved spectrum and the loss of potentially valuable structural informa- 
tion. Adding more coefficients gains nothing because the signal -to-noise ratio 
is of order one at the minimum. In addition, truncating the coefficients at 
the point where the white noise begins to dominate forces the dispersion of 
f(x) to remain finite. 

With this brief introduction it is now possible to list the six steps 
which together constitute the FFT/power spectrum deconvolution method: 

1 ) Both the raw data h(x) and the resolution function r(x) are 
background corrected. 

2) The functions h(x) and r(x) are Fourier transformed using an 
appropriate FFT algorithm. 

3) The Fourier coefficients of h(x) are examined and those up to and 
including the first minimum in the power spectrum of h are 
retained, the rest are replaced by zeros. 

4) The Fourier transform of f(x) is obtained by the point by point 
division of the Fourier transform of h(x) by the transform of r(x). 

5) f(x) 1s obtained by the inverse FFT. 

6) The solution is checked by forming the convolution of f(x) with 
r(x) and comparing the result with h(x). 

Each of these basic steps will now be described in further detail. 
Semiconductor detectors usually produce asymmetric Y-ray peaks with a non- 
linear underlying background. This is caused primarily by three processes: 
1) incomplete charge collection, 2) electron escape from the active volume, 
and 3) scattering of Y-rays (Chaglar et aU 1981). This background component 



41 

the deconvolution procedure is performed. The background correction routine 
employed in this work Involves a least squares fitting of a linear function in 
both the low and high energy tail regions. The region on the low energy side 
of the peak spans from 486 to 497 KeV, while the high energy region extends 
from 534 to 542 KeV. The straight line segments fitted in these regions are 
extrapolated to b03 and 518 KeV respectively where they are joined by a third 
line segment to produce a simple background approximation (Jackman et al . 
1974). The linear least squares analysis algorithm used in this work is a 
simplification of a higher order least squares routine by Gerald (1978). As 
in the work of Jackman et al. (1974), use is made of only the high energy half 
of the background corrected spectrum since it is known that incomplete charge 
collection of solid state detectors results in low energy tailing. In theory 
the spectrum should be symmetric about 511 KeV since statistically one should 
observe equal numbers of red and blue Doppler shifted Y-rays. 

Before proceeding to the details of the FFT algorithm one important 
point requires explanation. In order to correlate MCA channels with a cor- 
responding energy range, one must know the MCA energy resolution per channel 
as well as which channel corresponds to the center of the energy distribution 
and hence bll KeV. The former value is determined as described in the experi- 
mental procedures. However, the center of the distribution must be calculated 
for each experimental spectrum. A simple centroid calculation is unacceptable 
since the higher background on the low energy side would bias the result. 
Similarly t simply selecting the channel with the maximum number of counts is 
undesirable since, due tn the statistical nature of the experiment, the 
centroid channel and the peak channel are not necessarily the same. The tech- 
nique used in this work is to determine the eleven channel region which con- 
tains the lar ge st number of counts and call the middle channel of this region 
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the center of the distribution* In most cases, the channel corresponding to 
the center as calculated above was also the channel with the maximum number of 
counts. 

Proper arrangement of the data string 1s required before proceeding to 
the FFT algorithm. Since this algorithm requires periodic data, the raw data 
are rearranged in the following way* The first datum point of the rearranged 
set corresponds to the center of the raw data string. The remaining rear- 
ranged data points numbered from 2 to ((N/2) + l) correspond to the remaining 
values of the background corrected high energy half of the raw spectrum. It 
should be pointed out that N must satisfy the relationship N = 2" 1 where m is 
an integer. Since the real data string rarely has exactly ((N/2) + l) data 
values it Is usually necessary to insert a string of zeros at the end of the 
data set to fulfill this requirement. In order to achieve periodicity, data 
point N should correspond to data point 2. In general the points from ((N/2) 
+ 2) to N are defined in the following way: 

g((N/2) + P) = g((N/2) + 2 - P) (3.21) 
where P varies from 2 to (N/2). In the experiment performed in this work m 
was equal to eleven which corresponds to a total of 2048 data points (includ- 
ing the zerofill ). 

Once the functions h(x) and r(x) have been background corrected and 
properly arranged they are passed to the FFT algorithms. The FFT algorithms 
employed in this work are those developed by Cooley, Lewis, and Welch (1979). 
There are a number of other such algorithms available which should work 
equally well. In every case it is recommended that the Fourier coefficients 
be calculated with double precision accuracy. The advantage of using FFT's in 
place of discrete Fourier transforms (DFT) can be seen from the following 
5**R 1 ?? r .: -data- paints the--DFT T requireS-l£- oc-~ -4.2 x 10 6 . complex itbj1> . 
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triplications while the FFT requires 2N Log2 N or - 4.5 x 10 4 complex multipli- 
cations. In this example the number of calculations required is reduced by 
about two orders of magnitude. 

The justification for using only the Fourier coefficients up to and 
Including the one corresponding to the first minimum in the power spectrum of 
h(x) has been briefly discussed. Similar techniques are commonly employed in 
electronic signal processing and it can be shown that the number of coeffi- 
cients retained is a measure of the frequency resolution. In the present case 
the resolution refers to the width of the lineshape components which can be 
resolved. The procedure by which the computer calculates the number of re- 
tained Fourier coefficients from H(w) is as follows. Values associated with 
the power spectrum are obtained by multiplying each Fourier coefficient H(K) 
by its complex conjugate. However, since the function h(x) is even, the 
Fourier coefficients H(w) are all real. Truncating H(w) at the K correspond- 
ing to the first minimum in the power spectrum is equivalent to the require- 
ment that the set of retained coefficients be monotonical ly decreasing, which 
is easily programmed by testing 

H(w + 1) < H(w) (3.22) 
In this analysis only the first (N/2) + 1 Fourier coefficients are considered. 
The first coefficient corresponds to the D.C. component of zeroth frequency 
and the N/2 ♦ 1 coefficient corresponds to the highest frequency available. 
It is well known that the remaining coefficients correspond to the negative 
frequencies and are 1n fact the complex conjugates of the first (N/2 + 1) 
coefficients following the relationship g(N/2 ♦ P) = g*(N/2 + 2 -P) for P = 2 
to N/2. However since all the coefficients are real, no new information is to 
be gained by examining the negative frequencies. 
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Once the appropriate number of coefficients, say I, has been deter- 
mined the remaining N-I coefficients of H(w) are set equal to zero before 
dividing H(w) po1nt-by-po1nt by R(w) and proceeding to obtain f(x) by the in- 
verse Fourier transform. The consequences of this truncation, or filtering 
operation are well known. See, for example, the discussion by Verkerk (1982). 
In general the process results In some systematic errors in which the peak 
values are underestimated and nonphyslcal oscillations may appear In the de- 
convoluted function f(x). For Ooppl er-broadened positron annihilation spectra 
taken with our instrumentation this filtering error is negligible for Doppler 
spectra having more than 40K peak counts. 

To check the accuracy of the function f(x) requires the convolution of 
this function with r(x) and the comparison of the resulting function with 
h(x). The details of practical convolution techniques have been discussed by 
Lam etaU. (19»1), and the procedure used in this work is essentially the 
same. 

A complete listing of the Fortran FFT/power spectrum deconvolutlon 
program Is given in Appendix F. 

Test of the Method 

In order to explore the validity of the FFT/power spectrum method the 
following test was performed. The resolution function r(x) was obtained from 
the b!4 KeV gamma ray of 85 Sr, while the function h(x) was measured for a neu- 
tron Irradiated high purity single crystal sapphire sample. This sample had 
been irradiated in EBRII at 740°C to a fluence of 3 x 1<j21 neutrons/cm* with E 
< 0.1 MeV. Both r(x) and h(x) were then processed as described in the prev- 
ious section. For this test of the method f(x) was estimated from F(w) using 
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tral analysis. 

Figure 5 is a plot of the power spectrum for both F(w) and H(w). The 
filled circles represent the spectrum of the experimentally obtained function 
h(x). Note the rapid decrease of about six orders of magnitude in power with 
increasing frequency. This power spectrum is about what would be expected 
from a wide-band-width stochastic process contaminated with white noise. This 
white noise component is associated with the band of frequencies to the right 
of the first minimum in the power spectrum of h(x). Although the overall 
signal-to-noise ratio is very high, for frequencies in this white noise band 
the signal-to-noise ratio is less than one and therefore it is impossible to 
retrieve any structural information about the function f(x) in this high fre- 
quency range. This suggests that only those coefficients to the left of the 
white noise band, that is to the left of the first minimum, be retained. To 
use less than this number of coefficients in reconstruction is undesirable 
since it is known that the errors associated with truncation, or filtering, 
increase with a decrease in the number of coefficients used. 

The amplification of the statistical errors in h(x) by the deconvolu- 
tion procedure can be clearly seen in Figure 5. The open circles represent 
the spectrum of the deconvoluted function f(x). Note the drastic increase in 
the power associated with the high frequencies of f(x) as compared with the 
corresponding frequencies in h(x). However, the low frequency coefficients to 
the left of the first minimum are well behaved. These observations are the 
foundation on which the FFT/power spectrum method is based. 

It is not surprising to find that the number of retained Fourier re- 
construction coefficients, as determined by power spectrum analysis, is a 
function of the number of peak counts in the raw data h(x). Figure 6 is a 
semi -1 oy plot of this relatio nship f or i rradiated sapphire. This f i gure sug- 
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Figure 5. A Plot of the Power Spectrum of the Raw Data, h(x), (Filled 
Circles), vs. the Deconvoluted Data, f(x), (Open Circles). 
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gests that the relationship over this range of peak counts is an exponential 
of the form: 

y a m *n x (3.23) 
where y is the number of coefficients, x is the number of counts in the peak, 
and m is approximately equal to 7 for this data set. This curve demonstrates 
the high price one must pay in terms of data collection time to increase sig- 
nificantly the number of retained coefficients. Fortunately, our experience 
has shown that with the ~100 coefficients which are usually obtainable with 
40K counts in the peak (30 minutes of data collection) the results obtained 
are very satisfactory. 

Figure 7 is a comparison of the deconvoluted function f(x) and the raw 
data h(x) after background correction. The curves have been normalized to 
equal area. As expected, the deconvoluted curve is both sharper and smoother 
than the raw data. There appears to be a slight oscillation in f(x) In the 
tail region. The quality of this estimate of f(x) can be examined In the 
usual way by convolution with r(x) to obtain an approximation h'(x) to the 
original experimental function h(x) as is shown in Figure 8. The overall 
quality of the fit is so good that on this scale It is difficult to see any 
differences between the raw spectrum and the convoluted spectrum. The one 
exception occurs in the peak region where several of the raw data points of 
h(x) lie well below their counterparts in h'(x). Careful examination of the 
shape of h(x) in this region reveals some rather spurious changes in the slope 
of the spectrum. Physical constraints on the annihilation process suggest 
that, assuming no positronium formation, the annihilation spectrum should be 
relatively smooth in this region. This suggests that these slope changes in 
h(x) are related to statistical counting errors. To test this idea two other 
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One of the theoretical advantages of deconvolutlon is a significant 
improvement in the resolution of UBPAS measurements. Figure 9 is a plot of 
the high energy half of a typical 85 Sr peak before and after deconvolutlon. 
The closed circles represent the raw data after background correction, while 
the open circles show the shape of the resolution curve after deconvolution 
using 243 coefficients. This is the number of coefficients selected by the 
FFT/PS method corresponding to the first minimum 1n the power spectrum for 
this data set. Note that the full width half maximum (FWHM) value for the raw 
data is -1280 eV while the deconvoluted curve has a FWHM of -320 eV. This 
corresponds to a factor of four improvement in resolution. However, this 
result is somewhat misleading in that for aluminum the number of retained 
coefficients is approximately 60-70. The triangles in Figure 9 show the 
resolution function after deconvolution using 70 retained coefficients. In 
this case the FWHM is -576 eV corresponding to an improvement in resolution by 
a factor of -2.2. Clearly it can be seen that the improvement in resolution 
obtained by deconvolution using the FFT/PS method is a function of the number 
of retained coefficients, which in turn 1s directly related to the statistical 
quality of the raw data set (Schaffer et al. 1984). This Implies that the 
data collection process should be carried out for sufficiently long times in 
order to obtain good counting statistics. (The counting error is proportional 
to 1//N, where N is the number of counts in the channel). 

Unfortunately, counting for long times introduces another error 
source. This error source is related to instrumental electronic drift due, at 
least in part, to temperature fluctuations. There are two major ways in which 
other researchers have dealt with this problem. First, the accuracy of re- 
sults can be improved by operating the entire DBPAS system within a temper- 
ature CQnlrj^led_(^ l y^ used in 
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conjunction with the first, is to digitally stabilize the electronics. This 
digital stabilization is accomplished by simultaneously monltering the 511 KeV 
annihilation line along with a suitably located monoenergetic gamma ray source 
such as the 478 KeV gamma ray from he. In this procedure a spectrum is 
collected for -100S and then the centrold of the 478 KeV peak is calculated. 
Subsequently, several other spectra are collected and their corresponding 
centroids also calculated. Next, each spectrum is linearly shifted such that 
all of the 478 KeV centroids coincide. Finally the 511 KeV peaks are summed 
to produce a OBPAS spectrum with as many total integrated counts as is 
necessary to obtain suitable counting statistics. 

While temperature and/or digital stabilization should improve the re- 
producibility of UbPAS measurements, it is still possible to obtain meaningful 
results using unstabilized equipment. In this work neither of the stabiliza- 
tion techniques were employed, and the major consequence of such a procedure 
is an observable spread in the values of the lineshape parameters calculated 
for repeated OBPAS measurements of the same metallurgical sample. This prob- 
lem will be discussed in greater detail in the sections of this chapter deal- 
ing with the separation of multiple defect types. 

In order to minimize the change in the metallurgical defect structure 
of the samples as a function of time, all of the OBPAS measurements in this 
work were made at 77°K. Performing the measurements at liquid nitrogen 
temperatures raises an important theoretical question: is the assumption that 
the positrons are essentially at rest, with respect to the electrons, at the 
time of annihilation still valid? Previous experiments by Kublca and Stewart 
(1975) demonstrate that even at temperatures as low as 4.2°K the positrons 
reach thermal velocities prior to annihilation. 

sTZ7 - r fe-itetf Ffflnt-jJw • effect of-mea s u r -i ng -.■ the ann i h »a*4on „_ 
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lineshape at 77°K an experiment was performed in which six DBPAS spectra were 
obtained for a well annealed pair of aluminum samples* Three of the measure- 
ments were made at 77°K and three at room temperature* The small deviation 
(of order 1/VN) between the six spectra showed no correlation with measurement 
temperature. It was therefore concluded that decreasing the sample tempera- 
ture from room temperature to 77°K would not affect the positron-electron 
annihilation procedure in any experimentally observable manner* 

As mentioned in the experimental procedures the samples examined in 
this work were high purity aluminum disks heat treated to control defect type 
and concentration* It 1s assumed in the remainder of this work, without 
direct TEM evidence but based on the arguments presented 1n the experimental 
section, that one group of samples contained predominately monovacancies while 
a second group of samples possessed a defect structure dominated by prismatic 
dislocation loops. For each pair of disks three or four annihilation energy 
spectra were obtained, from which an array of lineshape parameters were sub- 
sequently extracted numerically. The lineshape parameters examined in this 
work include the first five even central moments of the energy distributions 
defined by the relationship: 




(4.1) 



where 



N 



is the number of channels in the spectrum 



x 



is the channel corresponding to the spectrum centroid 



f(x-j) is the number of counts in channel i 



and 



is an even integer (2,4,6,8,10) corresponding to the first 



through the fifth even central moments 



A second group of parameters is defined by the relationship: 
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i =x 

with j defined as j = I + 5(1-1) 

where f( Xi ) is the number of counts in channel i 

is the channel corresponding to the spectrum centroid 
and I is an integer from 1 through 12 

For example, x(l, is the ratio of the number of counts in the first five 
channel regions on the high energy side of the centroid (actually it includes 
the centroid) to the total number of counts in the hign energy half of the 
spectrum. Similarly. x(12) the ratio of the number of counts in the 
twelfth five channel interval (to the right of the centroid) to the total 
number of counts in the high energy half of the spectrum. 

The third, and final, group of lineshape parameters examined are the 
conventional "S" parameters (MacKenzie et al. 1970 ). These parameters are 

the total number of counts in the energy spectrum. The width of the energy 
window used to define the S parameter is selected so that the value of S is in 
the range 0.4 - 0.6 (Smedskjaer and Fluss 1983). m this study three 
S-parameters (S(1344). 3(1472), and S(1600)J were examined with respective 
central energy windows of 1.344, 1.472, and 1.600 KeV. 

These parameters were calculated for both the background corrected raw 
data and for the deconvoluted data. Hence a total of forty lineshape 
parameters were calculated for each OBPAS measurement. 

Figure 10 is a plot of the parameter x(l) for the deconvoluted data 
The open circles represent the values of this parameter for the "as quenched" 
samples (predominately monovacancies). The filled circles correspond to the 
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low temperature anneal samples (prismatic dislocation loops), and the tri- 
angles represent the samples after the high temperature anneal (the essen- 
tially defect free state). Note that the numerical value of x(l) increases 
with quenching temperature for both defect types. This trend is in agreement 
with theoretical predictions since an increasing defect concentration results 
in a higher probability of annihilation with lower kinetic energy conduction 
electrons, corresponding to a more narrow energy spectrum, and hence a larger 
value for x(l). The equations used to predict the theoretical defect concen- 
trations for both the vacancies and the prismatic dislocation loops as a func- 
tion of quenching temperature are discussed in Appendix 6. 

Figure 10 demonstrates that DBPAS measurements are sensitive to the 
presence of both vacancies and dislocations in high purity aluminum. This 
result is, of course, a necessary prerequisite for the possibility of DBPAS 
measurements being able to distinguish between the defect types. From this 
plot it is possible to perform some semiquantitative analysis and estimate 
some physical properties of aluminum. For example, one can estimate values of 
the "cross section" for the trapping of a positron by both vacancies and 
dislocations in aluminum. The accuracy of these estimates is limited by the 
reproducibility of the DBPAS measurements as well as the need for some 
independent information such as the lifetime of positrons as a function of 
defect type. 

In 1972 Cotterill et al_. published a paper in which they examined the 
average positron lifetime, T avg , on a set of aluminum samples very similar to 
the ones studied in this work. Their work clearly demonstrated that positron 
lifetime experiments were sensitive to both vacancies and dislocations in 
aluminum. Before discussing the similarities between their work and the 
current .experiments^ it j s_ important tp^ . 
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the parameters x avg and x(l). m a review article by Sieger (1980) it is 
shown that: 

T avg s (1 - At)r b + (4#3) 
where t 5 is the lifetime in the bulk material 
x t is the lifetime in the trapped state 
and A t is the probability that the positron annihilates from the 

trapped state 
Similarly, Sieger demonstrates that: 

F Ml - A t )F b ♦ AtFt (4 . 4) 
where F is any characteristic shape parameter obtained from the analysis 
of OBPAS measurements 
F b is the value of the parameter in the bulk state 
and F t is the value of the parameter in the trapped state 

Since the parameter x(l) can be shown to be equivalent to an "S" type 
parameter, and since an «S" type parameter is a common characteristic line- 
shape parameter. Eq. (4.4) can be rewritten as: 

X(l) « (1 - A t )x(l) 5 + Atxd^ (4 . 5) 
where x(l) 5 is the value of the parameter in the bulk state 
and x(l) t is the value of the parameter in the trapped state 
Clearly a comparison of Eq. (4.3) and Eq. (4.5) demonstrates the similarity 
between T avg and x(l). 

With this in mind it is not surprising to find that plots of t avg vs. 
quench temperature (T Q ) and x(l) vs. T Q share many common characteristics. 
(The x avg plot used for comparison is Figure 2 in the work of Cotterill et 
aK, 1972). In both figures the Tq at which the vacancy concentration is 
large enough to affect the measured parameter is ~200«C. The saturation value 
of both Javg *M*(l)Jor vwanc^occMJ? 
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the points corresponding to vacancies is sigmoidal in shape for both the 
lifetime and deconvolved DBPAS measurements. As for the points corresponding 
to samples with prismatic dislocation loops, the value of both parameters 
initially deviates from their bulk values for samples annealed after a quench 
from ~4UU°C. Neither parameter reaches a saturation value for dislocations 
within the temperature range examined In this work. 

Figure 10, therefore, shows not only the ability of deconvoluted DBPAS 
measurements to detect the trapping of positrons by both vacancies and dislo- 
cations in aluminum, but also qualitatively verifies the theoretically pre- 
dicted equivalence of the lifetime parameter r avg and the Doppler parameter 
x(D. 

Figure 11 1s a plot of the same llneshape parameter, x(l), vs. Tq for 
the background corrected raw data, i.e. the same data as shown in Fig. 10 but 
before the deconvolution procedure. Note the general similarity between these 
two figures. Careful examination of the raw data shows that the points 
corresponding to the 400°C quench are systematically displaced. That is, the 
entire group of five data points seems to be "too high" as compared with the 
rest of the data. To understand the cause of this shift it is necessary to 
examine the shape of the instrument resolution function R(x). 

During the course of these experiments the instrument resolution was 
periodically checked by collecting an energy spectrum for 8 5sr. The monoener- 
getic line is broadened due to the physical characteristics of the DBPAS sys- 
tem as previously described. This broadened energy line serves as the resolu- 
tion function, K(x), in the FFT/PS deconvolution procedure. It also describes 
the amount by which the undistorted energy spectrum is broadened by the in- 
strumentation. 
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Figure 12 shows two such resolution spectra. The open circles repre- 
sent the instrument resolution corresponding to the OBPAS measurements for the 
samples quenched from 400°C. The filled circles represent a resolution func- 
tion which is characteristic (within statistical error) of the resolution 
functions for all of the remaining data points. This comparatively large 
shape change in the function R(x) was apparently caused by a power failure in 
the laboratory. Note that the Function R(x) corresponding to the 400°C 
samples (henceforth termed SR850LD) is more sharply peaked, and therefore more 
narrow, than the other resolution function (SR85NEW). 

The physical significance of the change in resolution function is as 
follows. An observed energy spectrum will be broadened in two ways. The 
first way is by a decrease In the concentration of vacancy type defects, in 
the sample being examined. The second broadening mechanism 1s, of course, the 
broadening due to the instrumentation. Now, if the second source of broaden- 
ing remains essentially constant, then changes in the width of the observed 
spectra may be attributed to changes in defect concentration. On the other 
hand, if the broadening due to the instrumentation decreases, and this error 
source is not properly removed, the observed spectrum will be artificially 
narrow and its apparent defect concentration will be too great. Thus, if the 
data for the 400°C samples are collected with an artificially narrow resolu- 
tion function, i.e. SR8501D, their apparent defect concentration will be 
greater than they should be. 

Of particular significance is, that while a change in experimental 
conditions can alter the shape of the resolution function and thereby distort 
the values of lineshape parameters calculated for the raw data, the deconvo- 
lution routine, when supplied with the appropriate R(x) function, auto- 
^aMfJLliy _r!«l.d.s the desired_results 8 . _Thls, ilJusjtrat es that, for . reprpduca-,.. 
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bllity the deconvolution of DBPAS measurements should be considered. In 
addition, as mentioned previously, since deconvolution removes the instru- 
mental effects it should yield lineshape parameters which are directly 
comparable between different laboratories. 

Despite the advantages of deconvolution it 1s important to remember 
that this technique does, of course, have its limitations. In DBPAS 
measurements there are several sources of experimental error, each of which 
has a corresponding time scale. The magnitude of this time scale is the 
overriding factor in determining the applicability of deconvolution techniques 
in minimizing the specific error in question. For example, any change in 
experimental conditions which can be treated as a step function, that is, by 
changing from one constant value to another constant value in a very short 
period of time, can be easily treated by deconvolution. Such "error" sources 
include the accidental (e.g. power failure) or intentional changes in the 
electronic amplification of the energy signal. To correct for such a change 
one simply needs to collect a resolution function before and after the event 
and use the appropriate R(x) in the deconvolution algorithm. 

On the other hand, error sources having a time scale approximately 
equal to the time required to collect a data set (tens of minutes) can not be 
corrected by deconvolution alone. Examples of such possible error sources 
include line voltage fluctuations, temperature variations, and electronic 
drift. These error sources may cause the centroid of the resolution function 
to vary during the collection period, resulting in a composite R(x) which is 
artificially broad. In addition, even if the raw data set to be deconvoluted 
and the function R(x) are collected in rapid succession the amount of centroid 
drift for the two energy spectra may well be uncorrected. This type of error 
source requ1res_.a..djffe.cent : means of correction. The DBPAS measurements . 
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should be performed in a temperature controlled laboratory. Also, digital 
stabilization, as described previously, should be used to control centrold 
drift. 

The final error type is the small time scale or hiyh frequency class. 
An error source such as sixty-cycle electrical noise, if it were present, 
should average to zero over the time required to collect a data set. In 
addition, a high frequency error source should have a similar effect on both 
the resolution function and the DBPAS energy distribution. 

In order to judge the effectiveness of deconvolution in eliminating 
experimental error it is necessary to define some quantitative measure of this 
error. The predominant manifestation of experimental error appears as a 
spread in lineshape parameter values calculated for multiple DBPAS measure- 
ments on the same samples. Figure 13 is a plot of the parameter x(l) vs. Ty 
for the deconvoluted data, note the spread in values of x(l) for a given 
defect configuration. 

In an attempt to evaluate the performance of the deconvolution algor- 
ithm the following parameters are introduced (see Figure 13 for a physical 
interpretation). 

1) AP V /P 0 is a measure of the sensitivity of this parameter to the 
existence of vacancies. 

2) 6P V /AP V is a measure of the error associated with the spread in 
values of this parameter for vacancies. 

3) A Pi/ p o 1s a measure of the sensitivity of this parameter to the 
existence of dislocations. 

4) SPjyaPj^ is a measure of the error associated with the spread in 
values of this parameter for dislocations. 
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Clearly, it 1s desirable for the sensitivity (*P/P 0 ) to be as large as 
possible and the error (6P/AP) as small as possible for both defect types. A 
comparison of these quantities calculated from the deconvolved and the raw 
data is shown in Table 1. Note that for all four of the lineshape parameters 
listed the sensitivity ratio is higher and the error ratio is lower for the 
deconvoluted data than for the raw data. This analysis encourages the further 
use of the deconvolution procedure to improve the quality of the results ob- 
tained from OBPAS measurements. 

A tentative explanation for the apparent improvement in sensitivity 
and accuracy of the deconvoluted results is based upon trends in the shape of 
the energy spectrum before and after the application of the FFT/PS analysis. 
Figure 14 shows the high energy half of the peak region for an A£ sample 
quenched from ~60U°C. The deconvoluted spectrum is, of course, more narrow 
than the raw data but note that it is also smoother (due to the filtering of 
the high frequency components). One explanation for the comparatively large 
error ratios of the raw data is that they are related to the statistical 
counting errors in the peak region of the spectrum. If this is indeed the 
case it is not surprising to observe a decrease In the error ratio after the 
filtering of the high frequency components. However, the most probable 
explanation for the Increased sensitivity resulting from deconvolution is 
related to the removal of the "broadening" action of the experimental instru- 
mentation. Recall that the major effect of the instrumentation is to spread 
counts from one channel into several adjacent channels. Therefore, the true 
lineshape changes in the annihilation energy distribution, caused by the 
trapping of positrons at defects, are spread out over a much larger number of 
channels 1n the raw data. The deconvolution procedure corrects for this 
effect and in so doing. restores., some of the details of. the- "true" annihilation 
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energy dlstrlhutlon. The net result I. .„ Weas. ,„ defect sensitivity. It 

" e TO " ti0 " e '' ^"0" °> tne enh.nce.nent In sensitivity ana 

accuracy due to devolution noted In this work are consistent. Out sn-11. 
However, fro. this work alone It is not possible to estate the effect of 
t«nperature and electronic drift upon the results. !t would he e»tre»ely 

enlightening to repeat this analysis employing a temperature controlled and 

digitally stabilized system. 

Before moving on from the topic of parameter sensitivity and accuracy 
one more pertinent observation should be discussed. Figure 15 is a plot of 
the error ratio. *P /AP . calculated separately flt ^ 

(the parameter being examined is x(l». Several trends are apparent. First, 
It is clear that the error associated with the dislocations is larger than the 
error associated with vacancies. This trend may be a. direct consequence of 
the si 2e of the denominator. which is smaller, at a given quench tempera- 
ture, for samples with dislocations than for the corresponding samples with 
vacancy type defects. This observation is however somewhat deceptive in that 
the horizontal axis of Figure 15 is T Q and not defect concentration. A TOan - 
ingful comparison of parameter sensitivity between defect types should be made 
using equivalent defect densities. Unfortunately, the dislocation densities 
which can be introduced by this thermal treatment are not large enough for 
direct comparison. In addition, it should be mentioned that the trapping of 
positrons by dislocations may be somewhat more complicated than trapping by 
vacancies. This idea is discussed in detail by Martin and Paetsch (1 9 72). 
Very briefly, it may be that the depth of the potential well associated with a 
Prismatic dislocation loop is such that detrapping of the positron, due to 
local thermal excitations, becomes significant. If this is indeed the case, 
the interpretation of QBPAS data for, low, prismatic dislocation loop- densities 
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becomes more complicated. Clearly, additional experimental data 1s required 
to determine the exact reason for the observation that the error parameters 
are laryer for the samples containing primarily dislocation loops than for the 
samples dominated by vacancy type defects. 

Another interesting feature of Figure IS is the relatively constant 
value of 6P/AP for vacancies. This implies that the spread in the data de- 
creases with the quenching temperature, since clearly AP decreases with T Q . 
One possible explanation for this observation 1s related to the shift in the 
centrold of the experimentally obtained spectrum due to temperature and/or 
electronic drift. The explanation is that a centroid shift changes the peak 
shape of a broad and flat spectrum (low defect concentration) to a lesser ex- 
tent than that of a narrow and sharp spectrum (high defect concentration). 
This 1s. of course, an extremely qualitative argument but it does predict that 
in a temperature/digitally controlled series of experiments the relative size 
of the error bars associated with the high defect concentration samples should 
decrease. 

The preceding discussion has demonstrated the utility of deconvolution 
in DBPAS measurements. In summary, the use of the FFT/PS method improves the 
sensitivity and the relative error of the calculated Uneshape parameters. In 
addition, deconvolution automatically corrects for a wide range of error 
sources which result in changes in the instrument resolution function. Decon- 
voluted lineshape parameters should be more directly comparable between dif- 
ferent laboratories and, in addition, offer the potential ability to resolve 
fine structure in the energy spectrum which Is unavailable from direct obser- 
vation of the raw data. There is strong evidence (see for example Hood and 
Schultz, 1980) suggesting that the use of both temperature stabilization and 
digital stabilization, to control .centrold drift. ..shouj,d..s1gnlf±cantjy,.improve . 
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It is Interesting to note the conditions under which the undeconvo- 
luted raw data has yielded reasonable results in the past. The main require- 
ment is, of course, experimental conditions which yield a nearly constant 
resolution function, R(x). As one would expect, temperature and digital con- 
trol as well as performing the series of experiments as close together in time 
as possible, both tend to help yield reasonable results for the raw data. In 
addition, the use of "S" type parameters, which are essentially area integrals 
over large energy windows, will be affected only slightly by the "spreading" 
action of the instrumentation. The major advantages one gains by deconvolu- 
tion are: 1) the ability to compare more accurately results for a wide var- 
iety of experimental conditions and, 2) the possibility of observing a more 
detailed structure of the positron -electron annihilation energy spectrum. 

The Separation of Multiple Defect Types 

The previous discussion has demonstrated that deconvoluted OBPAS 
measurements are able to detect the presence of both vacancies and disloca- 
tions 1n high purity aluminum. The evidence for this statement comes from 
controlled laboratory experiments using samples with reasonably well charac- 
terized defect distributions. In order for OBPAS measurements to become use- 
ful as an NOT tool in industrial applications it is important for the tech- 
nique to be able to distinguish between multiple defect types. The reason for 
this requirement is that in a practical application one may require informa- 
tion about a particular defect type while being comparatively unconcerned with 
other defect types. For example, if one is interested in hydrogen embrittle- 
ment of high strength steels the "masking" of vacancies by nascent hydrogen 
may be of primary importance while slight variations in dislocation density 
may have relative!*" Tittie^i^Sificance. " '" • 
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Before discussing the results of this Investigation to determine the 
ability of DBPAS to separate defect types a brief review of previous work in 
this area will be presented. The basic concept is that no two defect types 
will have identical electron confiyurations associated with them. That is, 
the distribution of electrons (both in terms of electron density and momentum) 
in the neighborhood of a single vacancy will be different than in the neigh- 
borhood of a dislocation core. Qualitatively the atmosphere will be the same, 
i.e. a decrease In electron density and a lower average electron momentum, but 
quantitatively, the measured lifetimes and fine structure of the momentum dis- 
tribution should be different. 

This variation in electron distribution has been verified repeatedly 
in the case of lifetime measurements. For example, Fluss et aU (1978) report 
a lifetime of -244 ps for positrons trapped at a vacancy in aluminum while 
Hautojarvi et a]_. (1970) give the value of -229 ps for a dislocation in the 
same material. It should be mentioned that the accurate determination of 
multiple lifetimes from a single experiment requires some rather involved cal- 
culations. The complexity of the procedure increases rapidly with the number 
of unconstrained components one attempts to resolve and, for crystalline 
materials, the results are questionable if more than two defect types are 
considered. 

An alternate method of analyzing lifetime data is to calculate a 
single lifetime component, known as the "average" lifetime, for a sample con- 
taining multiple defect types. This approach simplifies the calculations and 
finds its greatest utility in experiments in which the concentration of one 
defect type is varying over a wide range of values while the concentrations of 
the other defect types remain relatively constant. Of course, such an 
— approach results in the loss of the abi 1 ity to separate defect types v 
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In the case of OBPAS measurements there have been fewer attempts to 
distinguish between different defect types. The most notable exception is the 
work of Mantl and Triftshauser (1978) which will be discussed in a following 
section. There are at least two major reasons for the lack of work in this 
area. The first reason Is related to the comparatively poor energy resolution 
of the OBPAS technique. Because of the smoothing action of the instrumenta- 
tion the fine structure of the energy distribution may have been hidden in 
previous experiments. In addition, even a slight change in the instrument 
resolution function during a series of experiments could mask the subtle 
changes 1n lineshape required to separate defect types. Although temperature 
and digital stabilization of the electronics would improve the situation, such 
a procedure does not address the issue of the smoothing action of the instru- 
mentation. It is the use of deconvoluted data which offers the possibility of 
examining the fine structure of the electron momentum distribution. 

As mentioned 1n the introduction, the most commonly reported Doppler 
lineshape parameter is the MacKenzie "S" parameter (1970). This lineshape 
parameter is defined as the ratio of counts In some central region of the 
DBPAS energy spectrum to the total number of counts in the spectrum. Since it 
is known that the energy spectrum becomes more sharply peaked with increasing 
vacancy type defect concentration it follows that the numerical value of "S" 
should be proportional to the total vacancy type defect concentration. The 
disadvantage of using this lineshape parameter is that it is unable to 
distinguish between multiple defect types. For example, if both monovacancies 
and edge dislocations are present in a sample both defect types will tend to 
increase the value of "S", however, it will be impossible to separate the 
relative defect concentrations solely from the calculated value of "S". 
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Mantl and Trifthauser (1978) developed a method for separating the 
effects of multiple defect types from experimentally obtained DBPAS data. 
To accomplish this they introduced a new lineshape parameter, I c , defined 
as the ratio of the sum of the number of counts in two regions, (one on the 
low energy side and one on the high energy side of the peak) divided by the 
total number of counts in the spectrum. This lineshape parameter, which in 
this work will be referred to as the winy or "W" parameter, can be shown to 
decrease with Increasing vacancy type defect concentration. When analyzed 
separately the "W parameter suffers from the same disadvantages as the "S" 
parameter in that it is incapable of separating multiple defect types. How- 
ever, Mantl and Trifthauser recognized that the combined use of "S" and "W" 
could yield useful results. 

In Appendix H the derivation of their defect specific M R" parameter is 
recreated. To summarize, it can be shown that in theory the "R" parameter, as 
defined below, is concentration independent, but characteristic of the type of 
trapping site Involved (Mantl and Trifthauser 1978). The "R" parameter is 
defined as: 

R - |(S - S B )/(W - W B )| (4.6) 
where S is the "S" parameter for a given specimen, 

S B is the value of "S" corresponding to a "defect free" specimen, 
U is the "W" parameter for a given specimen, 

W B is the value of "W" corresponding to a "defect free" specimen. 
While such a parameter is very useful, because of its ability to distinguish 
between multiple defect types, it still does not solve the problem completely 
due to some inherent theoretical limitations. The main restriction is that 
the analysis is valid if only one type of trapping center is present, or if 
one defect type is essentia) \\y _ dominant (Mantl aM_ Trifthauser 1978). In 
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practice this means the experimenter has the ability to calculate defect type 
and concentration provided each sample has only one predominant type of 
positron trap. The difficulties arise if a sample contains two defect types, 
say vacancies and dislocation cores, in roughly equal concentrations. 

In terms of the present experiments on high purity large grain size 
aluminum, because care has been taken to produce samples with only a single 
dominant defect type, one would expect the "R" parameter analysis to yield 
useful results. A family of K-type parameters can be calculated for the 
experiments in this work by examining ratios of the form 

RN = l((Ax(1)/Ax(j))| (4t7) 
where R N is an R-type parameter 

Ax(i) 1s the difference between x(i) for the current sample and the 
value of x(i) in a defect free sample 
and i is less than j 

In the original work of Mantl and Trifthauser (1978) the energy windows 
selected to define "S" and "W" were chosen so that the annihilation events 
contributing to the "S" parameter corresponded to conduction (valence) 
electrons while core electron events resulted in contributions to the "W" 
parameter. This relationship restricts the possible values of 1 and j as 
defined in Equation (4.7) above. 

The specific "R" parameter used in this study is defined as 

d _ I Ax(l) | 

1 " Kx(10) * x(ll) + x(12)j| (4.8) 
and is defined such that the energy windows are fairly close to those used by 
Mantl and Trifthauser (1978). For comparison. R x has a low energy window of 
-511 - 511.32 KeV and a high energy window of -513.88 - 514.84 KeV compared to 
windows of -511 - 511.5 KeV and -513.6 - 516.1 KeV in the referenced work. 
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Figure 16 1s a plot of R\ as a function of quenching temperature for 
the background corrected raw data. Note that although samples with primarily 
dislocation type defects have, in general, a higher Rj value than the samples 
with primarily vacancy type defects, there is a significant amount of scatter 
in the data. The source of this scatter is not known exactly but it is prob- 
ably related to the fact that the data was collected without temperature or 
digital stabilization. In fact, it is plausible that the spread in Ri values 
is a numerical amplification of the spread in values of x(l), x(10), x(ll) and 
x(12). In addition to this error source, any error associated with changes in 
the instrument resolution function, (IRF), should also be present in Figure 
16. To test the utility of deconvolution in R-parameter analysis one would 
like to have a data set in which the second error source is present but the 
first error source, related to electronic drift, is negligible. Unfortunate- 
ly, such a data set can not be obtained with the instrumentation used in this 
study. 

In an attempt to isolate the instrument resolution effects the fol- 
lowing procedure was employed. Using the fact that all DBPAS measurements for 
a given metallurgical sample were performed in a comparatively short time 
interval (~8 hrs.), it seems reasonable to assume that the spread in calcu- 
lated lineshape values is a result of the lack of temperature and digital 
stabilization rather than of changes in instrument resolution. This hypoth- 
esis is supported by the observation that for a given set of consecutive 
experiments, the IRF's measured before and after the experiments showed only 
very slight changes. If this assumption is true, one may be able to arti- 
ficially select a subset of the original data which qualitatively corresponds 
to a series of experiments with reduced electronic drift error. In this work 
each group of DBPAS measurements, corresponding to a sample subjected to _ .... 
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a given thermal treatment, was Independently examined and an average value for 
each relevant llneshape parameter was calculated. The corresponding Rj values 
calculated from the average llneshape parameter values are shown in Figure 17. 
These values are. as In Figure 16. for the background corrected raw data. 

A similar procedure was then employed for the deconvoluted data. That 
is. average x(i) values were calculated for each specimen and then used to 
calculate R x values. These R x values for the deconvoluted data are shown in 
Figure 18. It is quite apparent that the spread in R-parameter values has 
been increased by tne FFT/PS deconvolution algorithm. The explanation for the 
increased scatter in the data Is not immediately obvious. In a previous 
section of this work It was demonstrated that deconvolution improved the 
quality of the calculated x(l) lineshape parameter values. That is. the 
spread in x(l) values for DBPAS measurements performed on a given sample was 
decreased (refer to Table 1). One could argue that decreasing the scatter of 
x(l) for a given sample should improve, or at least not adversely affect, the 
amount of scatter in the calculated Rj values. The difficulty is, of course, 
that Ri depends not only on x(l) but also on the values of the parameters 
x(lO), x(ll), and x(12). These parameters are defined in the tail regions of 
the energy distribution and it is in this region that one must look to under- 
stand the difficulty associated with the calculation of R-type parameters for 
FFT/PS deconvoluted data. 

Figures 19 and 20 are a comparison of the energy distribution (raw 
data) for a sample quenched from 600°C and then subsequently annealed at 
523°K. Figure 19 shows the differences in the peak region of the spectra. 
Recall that x(l) is defined over the energy window of 511 - 511.32 KeV. 
In Figure 20 the vertical scale has been expanded so that the differences 
in the tail, region are visible. Several important, features.. of _tt»e_eneegy_ 
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distribution In the tan region are apparent. First, it must be noted that 1n 
general the data points for the quenched sample He below their counterparts 
for the annealed sample as predicted by theory. Secondly, the statistical, or 
counting, error from channel to channel can be recognized. The former obser- 
vation 1s believed to account for the qualitative agreement between the- 
oretical predictions and the calculated Rj values for the raw data. The 
latter observation, statistical noise, is the most probable explanation {along 
with electronic noise) for the spread in Rj values as shown in Figure 16. 

For comparison. Figures 21 and 22 show the change in energy distribu- 
tion for the same two samples after deconvolutlon using the FFT/PS algorithm. 
In Figure 21, which shows the peak region, some of the desirable features of 
deconvolutlon are evident. Most notably, the statistical channel to channel 
noise has been removed (compare with Figure 19). A careful comparison of the 
raw and deconvoluted data shows that deconvolutlon does indeed decrease the 
full width half maximum as discussed previously. 

Figure 22 shows, with an expanded vertical scale, an undesirable 
feature of the FFT/PS deconvolutlon method. It was mentioned in Chapter III 
that truncating the Fourier coefficients is equivalent to the insertion of a 
low pass filter and it 1s well known that low pass filtering Introduces 
certain errors. Recalling that the Fourier analysis of the convolution 
equation, neglecting the statistical counting error, yields the following 
equation for the desired function f(x): 



f(x) = / u(w)F(w)e- iw * dw (4ao) 



where u(w) = 1 for w < w c 
= 0 for w > w c 
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and w c corresponds to the frequency which is the first minimum of the 

power spectrum of h(x) 
It can be shown (Bendat and Piersol 1971) that error introduced by such a 
filter is proportional to (l/w c )sin(2wxw c ). The significance of this result 
is that the amplitude of the error is inversely proportional to w c (the number 
of retained coefficients), and the frequency of the resulting error is 
proportional to w c . 

This effect is clearly evident in Figure 22. The open circles 
represent the annealed sample for which the power spectrum analysis yielded a 
w c of fifty six (56) coefficients. For comparison, the filled circles, 
corresponding to the quenched specimen are calculated using sixty seven (67) 
retained coefficients. It is quite obvious that the error amplitude is 
decreased for the sample with the larger number of retained coefficients. A 
careful analysis of the data shows that the frequency of the error is 
proportional to w c . 

With this information, it is now possible to explain the increased 
scatter in Rj values displayed by the FFT/PS deconvoluted data. Because the 
energy window used to define the wing parameter (513.88 - 514.84 KeV) is in 
the range where the amplitude of the filtering error is of the same order of 
magnitude as the amplitude of the true signal, this error source is 
significant. If the wavelength of the error source was small compared to the 
size of the energy window, then the error contribution would average to nearly 
zero and this difficulty would not exist. The problem arises because the 
eneryy window spans less than a period of the sinusoidal error source and it 
spans a different fraction of the period depending on the number of retained 
Fourier coefficients. For example, in Figure 22 the energy window of interest 
for the annealed sample spans approximately the, positive half cycle of _the 
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sinusoidal filtering error. This leads to a calculated wing parameter with an 
artificially high value, which In turn distorts the value of the calculated Rj 
parameter. Since multiple 08PAS measurements on even the same metallurgical 
sample may yield slightly different numbers of retained coefficients, this 
error source can lead to a scatter of Rj parameters for a single sample. 
Because the number of retained coefficients varies more widely between samples 
of differing defect density (or defect structure), the spread in Rj values 
becomes larger for such samples. 

Having related the scatter In Ri values for the deconvoluted data to 
the phase and amplitude of the filtering error, it is now possible to try to 
correct for this error source. There are at least three ways to proceed. 
First, one could simply extend the experimental collection times to get a 
higher number of total integrated counts In the energy spectrum. This would 
result in better counting statistics, and a higher number of retained 
coefficients (Schaf fer et aK 1984). Since an Increased number of retained 
coefficients decreases the amplitude of the filtering error, such a procedure 
should yield better Ri values. Unfortunately one must pay a very high price 
in terms of collection time to significantly increase w c . 

Another approach is to use the same number of retained coefficients 
for each data set. That is, one could calculate a w c value for each data set 
to be analyzed and then deconvolute all of the data sets using the minimum 
value of w c . Using this procedure will not eliminate the error source, but it 
will be controlled so that Its contribution to the calculated values of Rj 
does not vary between experiments. Clearly such an approach has several 
inherent disadvantages. The primary disadvantage is that the absolute values 
of Rj will be a function of the particular value of w c used in the deconvolu- 
tion urocedure. One would prefer, an approach which eliminates, or at least 
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greatly reduces, the filtering error rather than simply controlling it. A 
method for accomplishing this goal will be discussed in the final chapter of 
this work. What follows here is a discussion of controlling the error source 
1n the tail region by using the same number of retained coefficients for each 
data set. 

Figure 23 shows the energy window used to calculate the wing 
parameter. The vertical scale has been greatly expanded to show the fine 
detail of the structure. This plot is for the same samples as in Figure 22, 
but in this case both data sets have been deconvoluted using fifty seven (57) 
coefficients. As expected, the filtering error is in phase. In addition, as 
theory predicts, the data points corresponding to the quenched sample are 
below the points representing the annealed sample. One can extend this 
procedure to all of the other data sets and then calculate the corresponding 
Rl parameters in the manner previously described. The results are presented 
in Figure 24. When this figure Is compared with Figure 18 it can be seen that 
the filtering error does influence the scatter in R x values. By simply 
controlling this error source. I.e. by using the same number of retained 
coefficients in the deconvolution of all data sets, the quality of the 
calculated Rj values can be improved. 

To summarize this section, the average values of Rx» along with the 
associated standard deviations for each of the two defect types, have been 
calculated for four different types of data. The four data types are: 1) the 
entire raw data set (Fig. 16), 2) the average raw data (Fig. 17), 3) the 
averages of the deconvoluted data as obtained directly from the FFT/PS 
algorithm (Fig. 18), and 4) the averages of the deconvoluted data using the 
same number of retained coefficients "for each data set (Fig. 24). These 
values are shown in Table II. - 
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Table II 



Data Type 


Ravg< Vac -) 


Ravg(° is -) 


o(Vac.) 


o(D1s.) 


1 


.261 


.458 


.0677 


.1751 


2 


.240 


.440 


.0347 


.0492 


3 


.234 


.364 


.1515 


.1876 


4 


.253 


.424 


.0216 


.0290 



Note that the largest standard deviations occur for the deconvoluted 
data when the maximum number of coefficients are used in the reconstruction of 
each f(x). This result is explainable in view of the previously discussed re- 
lationship between w c and the amplitude and phase of the filtering error. 
When this error source is controlled by using a single value of w c for the de- 
convolution of all data sets (data type 4) the calculated standard deviations 
are improved. With this modification to the FFT/PS deconvolution algorithm it 
appears that if one defect type dominates, one can distinguish between 
vacancies and prismatic dislocation loops* in high purity large grain si2e 
aluminum. 

To conclude the results and discussion chapter of this work the even 
central moments of the annihilation energy distributions will be discussed. 
As mentioned in the introduction, because this technique has the potential to 
distinguish between multiple defect types, this family of lineshape parameters 
is interesting from a theoretical point of view. Briefly, this possibility 
arises because the higher order central moments are increasingly dominated by 
the tail regions of the energy distribution. Thus, the low moments give 
information about the central region (S-type information) while the higher 
order moments yield W-type information. Unfortunately, this analysis is 
extremely sensitive, to small .errors in the. tail regions which are greatly 
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amplified by moment analysis. Therefore the statistical error ruled out 
moment analysis in the raw data while filtering error precluded its use in tl 
deconvoluted data. It is possible, however, that an improved filtering 
technique, discussed in the following chapter, will allow moment analysis to 
become a useful tool in future work. 



CHAPTER V 
CONCLUSIONS AND FUTURE WORK 



Over the past several decades positron annihilation spectroscopy has 
evolved Into an accepted tool for the study of the defect structure of mater- 
ials. Of the three most common measurement techniques Doppler broadening of- 
fers the advantages of reasonably high data collection rates, comparatively 
low source strength requirements, and few geometrical experimental con- 
straints. It is therefore the technique with the best chance of making the 
transition from a laboratory research device to an industrial NDT tool. Un- 
fortunately. DBPAS measurements suffer from some inherent limitations. The 
two most significant limitations are comparatively poor instrumental energy 
resolution and limited ability to distinguish between multiple defect types. 
Based on these observations, the two primary goals of this work were: 1) to 
develop an efficient, model independent, noniterative deconvolution method, 
and 2) to improve the ability of DBPAS to distinguish between multiple defect 
types. 

In order to improve the effective energy resolution of the DBPAS in- 
strumentation the fast Fourier transform/power spectrum analysis deconvolution 
algorithm was developed ( Sc ha f fer et al_. 1984). The convolution integral 
equation was shown in Chapter III to describe the relationship between the ex- 
perimentally observed energy distribution, the instrument resolution function, 
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data is that results obtained in different laboratories should be more direct- 
ly comparable if the instrumental broadening effects have been accurately re- 
moved. This ability to directly compare numerical values for standard line- 
shape parameters should help to clear up some of the discrepancies which 
currently exist in the literature. 

While it is fairly clear that accurate deconvolution should correct 
for changes in the IRF and improve the agreement of results obtained using 
different instrumentation, some of the other trends observed from the compar- 
ison of lineshape parameter calculated before and after deconvolution were not 
entirely predictable. For example, based on the limited results presented in 
this work, both the sensitivity and the accuracy of the lineshape parameters 
are improved by deconvolution. One explanation for the improved accuracy, 
i.e. the decrease in the error ratio, of the deconvoluted lineshape parameters 
1s that the filtering operation has essentially removed the statistical count- 
ing errors present in the raw data. A comparison of multiple DBPAS measure- 
ments performed on the same sample show that before deconvolution, differences 
in the statistical counting error 1n the peak region can lead to a spread in 
calculated x(l) values. However, after deconvolution, the energy distribu- 
tions show less variation from experiment to experiments, and hence yield x(l) 
values with correspondingly less scatter. 

The increased sensitivity of the deconvoluted parameters may be 
related to the removal of the "broadening 1 ' action of the experimental instru- 
mentation. Recall that the major effect of the instrumentation is to spread 
counts from one channel into several adjacent channels. Therefore, the true 
lineshape changes in the annihilation energy distribution, which are caused 
by the trapping of positrons at defects, are spread out over a much larger 
number of channels in the—raw data* The deconvolution .pxQ<Le^re.^GOEte^Lt3^f oc 
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this effect and, In so doing, restores some of the details of the "true- 
annihilation energy distribution. This results in an increase in observed 
defect sensitivity. 

Even with the improvements in lineshape sensitivity and accuracy due 
to the FFT/PS algorithm, the DBPAS technique must have the additional ability 
to distinguish between multiple defect types if it is to become a useful tool 
in the study of industrial materials. In order to investigate the usefulness 
of deconvolutlon in resolving multiple defect types, a comparison was made 
between the values of the R x parameters calculated before and after deconvolu- 
tlon. If the number of reconstruction Fourier coefficients is allowed to vary 
between samples, as would be the case if the FFT/PS method is used in its 
standard form, then the calculated R x parameter value yields little useful 
information. 

In order to explain this somewhat negative result, the consequences of 
truncating the Fourier coefficient, also known as low pass filtering, were ex- 
amined in greater detail. It is known that the truncation of the coefficients 
introduces a filtering error with an amplitude Inversely proportional to the 
number of retained coefficients, w c , and a frequency which is proportional to 
w c . Ideally, one would like to greatly increase the number of retained coef- 
ficients in order to reduce the amplitude of this filtering error. Unfortu- 
nately, as w c is increased beyond the first minimum of the power spectrum the 
calculated solution will contain large, rapidly oscillating components which 
are not physically plausible. The mathematics of this argument were presented 
in detail in Chapter III. 

Since it is not possible to reduce the amplitude of the filtering 
error by arbitrarily increasing w c an alternative method was developed to 
improve the calculated values ofLRj fordeconvoluted data. The first step in- 
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the procedure Is to calculate a w c value for each of the data sets to be 
analyzed. Next the minimum value of w c is selected and subsequently employed 
in the deconvolution of all the data sets. Using this procedure does not 
eliminate the filtering error, but rather controls the source such that its 
contribution to the calculated R x values remains approximately constant from 
experiment to experiment. Although this solution is not ideal, the data 
presented in Chapter IV clearly indicate that by controlling the error source 
the reproducibility of the R x parameter, and Its usefulness in the separation 
of multiple defect types, is improved. 

To summarize this section it can be stated that the FFT/PS deconvolu- 
tion algorithm appears to improve the quality of the data as obtained from 
UBPAS measurements in several ways. These improvements, as presented In this 
work include: 1) automatic correction for accidental and/or intentional 
changes in the instrument resolution function, 2) the potential of direct 
comparison of quantitative results between laboratories, 3) an improvement in 
both lineshape sensitivity and accuracy, and 4) the ability to distinguish 
between vacancy and dislocation type defects, as demonstrated for large 
grained high purity aluminum. 

To generalize these conclusions, further studies of this kind should 
be performed using a variety of different alloys with known defect type and 
concentration. In addition, it would be extremely useful to repeat the entire 
series of experiments using a DBPAS system which is both digitally and temper- 
ature stabilized. Such a data set would allow one to separate the effects of 
instrument broadening, which are theoretically correctable by deconvolution, 
from the nonsystematic effects of instrumental drift. In addition, spectra 
collected on a stabilized system could conceivably possess better counting 
statistics (due to longer allowable experimental t InftsjL-JitLl eh could in . turn ^ ■ 
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increase the number of coefficients retained in the deconvolution procedure. 

Additional experiments should be performed on the FFT/PS algorithm. 
These experiments could be designed to investigate alternative filtering 
schemes with the hope of reducing the associated filtering error. Most 
notably, the use of cosine tapering appears to be a viable approach. Once the 
filtering error has been reduced or eliminated the utility of moment analysis 
in the separation of multiple defect types should be reinvestigated. It 1s 
not known at this time if the current FFT/PS algorithm is capable of disting- 
uishing the formation of positronium in a nonmetallic material. Pos1tron1um 
would show itself as a relatively high frequency characteristic of the annihi- 
lation energy distribution. The key factor is determining whether the fre- 
quencies required to observe positronium are to the left of the first minimum 
in the power spectrum. This question could be investigated using both numeri- 
cally simulated data and data obtained from materials in which positronium Is 
known to form. 

In order for Ooppler broadened positron annihilation spectroscopy to 
become a widely used industrial NOT tool it must continue to evolve in several 
ways. Its energy resolution must improve, the ability to separate multiple 
defect types must be further developed, and the reproducibility and reliabil- 
ity of the experimental results must increase. Althouyh the fast Fourier 
transform/power spectrum analysis deconvolution technique has clearly not 
solved all of these problems, the results presented in this work indicate that 
this method measurably improves the quality and reproducibility of the infor- 
mation obtainable from DBPAS measurements, and serves as a sound basis for 
future refinements. 
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APPENOIX A 



In this appendix the equations used to theoretically determine the 
stopping* or thermal izati on, times of positrons in metals will be discussed. 
For a more complete treatment of this subject see Brandt and Arista (1982). As 
briefly mentioned in the text, positrons slow down mainly by two mechanisms: 
1) ionization and excitation of electrons which is characterized by the elec- 
tronic stopping power, S e » of the material, and 2) phonon interactions charac- 
terized by the phonon stopping power, S p ^. The second mechanism dominates at 
epithermal energies below some critical value E c characteristic of the mater- 
ial. In the range from the maximum source energy, E s , to E c the electronic 
stopping power dominates. This energy range can be further divided Into two 
subranges separated by the photon energy for which the material has its maximum 
stopping power, E m . In the range E |n < E < E s an accurate estimate for the 
stopping times is given by 

t! » 7.7p-l(E i /mc 2 ) 1 - 2 (A.l) 

where t] is the stopping time in psec 

p Is the density of the material in g/cm 3 
mc 2 is O.bll MeV 
and Ef is the initial positron energy 
since 

Ej/mc 2 < E s /mc 2 - 1 (A. 2) 

we can rewrite Eq. (A.l) as 
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*' " 7 - 7 °-' (A.3) 
-hich for „ st „ t , )s of jMerest resu]ts u Mopp)n9 tiMs ^ _ ) ^ ^ 

high energy region. 

For E c < £ < £„, the electrow gas theopy Qf conduction electrons w . th 

Fermi energy E F is used to produce the equation 

S = -Lr "I f2i3/2 E 5/2 

e 10b h l-m 1 "Tf (A.4) 
in the case of a continuous slow down approximation 

t = - lfhl/2 r E f dS 
s l 2 J J E i JUTS (A.5) 

Substituting E«. (A.4) into Eq. (A.5) and integrating yields 

t .105 h£ F 
2 " Bit TT (A.6) 
L f 

where t 2 is In picoseconds. 
In the low energy region 3/2 KT < E < E c phonon interactions dominate with 

S =1 m 3 w 2 p 
PH « 77 E (A.7) 

where w is the lattice deformation potential. 
Eq. (A.4) can be combined with Eq. (A.7) to define E c yielding 

£ 3/2 _ 105b 2 m 5 / 2 E F 3 
° 2 ,/2 n2 ~~{?~~ (A.8) 

Inserting Eq. (A.7) fBto Eq . (A . 5) and lntegrat1 „ g from 3/2 KJ ^ ^ 

by Eq. (A.8) yields 
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using typical values for common metals it can be shown that 

*3 > X Z > *1 (A. 10) 

Thus, phonon thermal ization times are larger than electronic stopping times. 
The total stopping time is given by 

t $ - t! + t 2 ♦ t 3 (A. 11) 

with typical theoretical values of 13 psec in At and 45 psec in Pb (Brandt and 
Arista 1982). 



APPENDIX B 



The purpose of this appendix is to describe the relationships between 
the momentum distribution of the electrons in the sample and the measured 
Doppler shifts in energy as well as the angular deviations from colinearlty of 
the Y-rays. In the center of mass rest reference frame the annihilation T-rays 
will be emitted in directions separated by 180° and each with an energy of bll 
KeV which is equal to the rest mass of an electron (or positron). However, in 
the laboratory frame the momentum of the electron (the positron is assumed to 
be at rest as described in the text) produces changes in both the observed y- 
ray energy and the angle between the Y-rays. The geometry of the situation is 
shown in the diagram below, 

^ Center of Mass 
\ Velocity V 




Straightforward application of trie Lorentz transformations of special 
relativity combined with some trigonometric manipulation yield the result that 

«♦ - sin * Q + e (^) (8.1 ) 

where 6^ is the anjjular deviation from 180° 

v is the electron velocity 

(107) 
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c is the speed of light 

♦o is the angle between the electron velocity vector and the 
direction of y-ray emission 

y 2 2 
and e^j is a term of order which is negligible for v « c 

(For a detailed derivation of this result see Rindler 1960). The relativistic 
definition of momentum 1s 

P = am 0 V (B.2) 

where p is the momentum 

a = (1 - V)" 1/2 
c 

and m,, is the particle mass 
For v « c, o > 1 and p - mv 
Substituting this result into Eq. (B.l) yields 

6* - §B-. sin * o (B 3) 

Defining the transverse momentum component as 

PT * P sin (B A) 

and substituting Eq. (B.4) into Eq. (B.3) yields the final result 

2P T 

6 * "mT (B.5) 
o 

demonstrates that measuring the distribution of the angular deviations. 6*. is 
in essence equivalent to measuring the distribution of the momentum of the 
electrons within the sample. 

In the case of Doppler measurements one has 

6P Y = P e cos <» 0 (8>6) 
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where 6P Y 1s the change in momentum of the Y-ray 

P e 1s the electron momentum 
and <j> 0 is as defined above 
Using Eq. (B.2) with v « c Eq. (B.6) takes the form 

6P T = m^ cos * 0 ( B . 7 ) 
It can be shown that 1n relativistic dynamics 



(B.8) 



E 2 = P 2 C 2 ♦ ^2 
where E is the relativistic energy 
p is the relativistic momentum 
c is the speed of light 
and nio Is the particle rest mass 

Now for a photon (Y-ray) m,, = 0 so Eq. (B.8) reduces to 

E Y s Py c (B.9) 
Combining Eqs. (B.6), (B.7), and (B.9) yields 

«E Y * 6 Pt c = p e (cos * 0 )c - mgCVetcos ♦„) (B.10) 
where 6E T is the Doppler energy shift of the observed Y-ray 
6p Y is the change in momentum of the Y-ray 
and M e is the momentum of the annihilating electron 

From Eq. (B.lO) it can be seen that the momentum distribution of the annihi- 
lating electrons can be deduced from the measured Doppler energy shifts. 

Combining Eq. (8.5) with (B.lO) yields the following relationship 
between 6$ and 6E: 

6E = 255.56* (B>n) 

where 6E is in eV 
and <5$ is in mrad 
Finally, by combining Eq. (B.lO) and the fact that for an electron m 0 c2 = 511 
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KeV it can be shown that 

6E Y = 511 (cos * 0 )(-^) (B.12) 

where 6E Y is in KeV 

For the case of a 1 eV electron (v e = 5.93 x 10 7 cm/sec) moving in the direc- 
tion of Y-ray emission, Eq. B.12 predicts a 1.01 KeV shift in observed energy. 



APPENDIX C 



The material examined 1n this work was high purity (99.9995%) aluml 
supplied by the Materials Research Corporation. The chemical analysis, In 
terms of PPM of metallic interstitlals (and some gases) 1s listed below. 



Element 


PPM 


Element 


PPM 


El piripnt* 


PPM 

r rn 


Ag 


<0.1 


LI 


<0.1 


u 




Al 


major 


Mg 


1.5 


v 




As 


<0.1 


Mn 


<0.1 


u 


<\J 9 c 


Au 


<0.1 


Mo 


<0.1 


Y 




B 




Na 


0.2 


Zn 


0 4 


Ba 




MK 


<U.l 


Zr 


<0.1 


Be 




N1 


0.1 






Bi 


<0.1 


Pb 


<0.1 


CI 


0.6 


Ca 


<0.1 


Pd 


<0.1 


F 




Cd 


<0.1 


Pt 


<0.14 


H 


3 


Co 


<0.1 


Hh 


<0.1 


N 


5 


Cr 


<0.1 


Sb 


<0.1 


0 


20 


Cu 


0.1 


SI 


0.9 


P 


0.3 


Fc 


<0.1 


Sn 


<0.3 


S 




Ga 


<0.1 


Ta 








Ge 


<0.1 


Th 








In 


<0.1 


Ti 


<0.1 






K 


<0.1 


Tl 









(in) 



APPENOIX D 



In this Appendix the utility of Fourier analysis in the solution of the 
convolution integral equation will be demonstrated. In the first section of 
Chapter III it was shown that 

h(y) = / f(x)r(y - x)dx ( D .i) 
By definition of the Fourier transform 

F(w) = /" eiwx f(x)dx (D 2) 



Therefore 



F(w)K(w) = [/_ e1wx f(x)dx][ ^ e iwt r(t)dt] (Q 3) 

OO DO 

" J m i m e iwx f(x)e 1wt r(t)dxdt ( D .4) 

<B OO 

" L L eiW(x+t)f ( x ) r (Odxdt (0.5) 



Now let x + t = y 
so t - y - x 
and dt = dy 
Hence Eq. (0.5) can be rewritten as 



F(w)R(w) = / / e iwy f(x)r(y . x)dxdy ( 

OD OO 

= / e1wy C / f( x )r(y - x)dx]dy (D . 7) 
(112) 
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But the Integral within the brackets is by Eq. (D.l) equal to h(y). Therefore 
Eq. (U.7) can be rewritten as 

F(w)R(w) * / ei*yh(y)dy (0.8) 

and using Eq. (D.2) on the right hand side 

F(w)R(w) = H(w) (l). 9) 
This demonstrates that in Fourier space convolution is simply a complex multi- 
plication. 



APPENDIX E 



The purpose of this Appendix is to demonstrate some of the mathematical 
difficulties associated with the incorrectly posed convolution integral 
equation. The primary source of the problem is that the left hand side of the 
equation 

jt b 

* (x) = / f(y)r(x - y)dy (E.I) 
a 

is only known to within some error e. That is, if h(x) is the "true" function 
then 

max |h*(x) - h(x)| < e (E.2) 
a<x<b 

Any function fj (y) which when convoluted with the known resolution function 
r(x) yields a function hj(x) that satisfies Eq. (E.2) must be considered a 
solution of Eq. (E.l). 

Assume there are two solutions to Eq. (E.l) which are f](x) and f£(x) = 
fj(x) + c s1n(2imw/(b-a). By definition 

* b 

hi(x) - / f](y)r(x - y)dy (E.3) 
a 

and 

* b 

h 2 (x) 5 / if](y) + C sin (2wny/(b - a))]r(x - y)dy (E.4) 
a 

Subtracting Eq. (E.3) from Eq. (E.4) yields 
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* * b 

h 2 (x) - hjfx) = / r(x - y)C sin(2*ny/(b - a))dy (E.5) 



a 

Now integrating by parts with 



yields 



u ■ r(x - y) dv = C sin (|l^-)dy 

au . -(r(» y - W) ay , . . C£|^ ' os ^ 

/J r(x - y)C st„(|^j)d y - r(, - y) C- C '|,- cos (|^)^ 

Now integratlny the last term on the right hand side by parts with 

- - cos <fH*j) a> . "(r(« d; r)) dy 
"« ' " if 2 ?*!) =i» (f^) * « ■ r(x - y> 



Yields 



£ ~£?tJ "'^oy' "" ■» ' tr(x - y) cos 

- % - " yXf^) *<n <f^> dy (E . 7) 

Combining Eqs. (E.6) and (E.7) yields 

jj f (x - y) C sin (|^> dx = ^ ») E r(« - y) cos (£.8) 

and combining Eqs. (E.8) and (E.S) yields the result 

l£(x) - h^x) -| [(^) r(x - y) cos (f=?^)3j (E.9) 



So it can be seen that for any value of C one can select a value of N large 
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* * 

enough so that the difference between h 2 (x) and t\ } (x) is smaller than any 
previously determined number. Therefore, both h*(x) and h 2 (x) satisfy Eq. 
(E.Z) and f!(x) as well as f 2 (x) must be solutions to Eq. (E.l). The important 
point Is that one may select C and N in such a way as to allow f](x) and f 2 (x) 
to differ from each other by as much as one likes. Clearly the problem as 
posed does not have a unique solution and a different solution technique is 
required. 



APPENDIX F 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

C THIS IS A LISTING OP THE FORTRAN 77 VERSION OF THE PROGRAM C 

C DECONFFT WHICH IS DESIGNED TO DE CONVOLUTE DOPPLER BROADENED C 

C POSITRON ANNIHILATION SPECTRA. FOR A DESCRIPTION OF THE C 

C MATHEMATICS INVOLVED SEE J. P. SCHAFFER, E.J. SHAUGHNESSY, C 

C AND P.L. JONES, JOURNAL OF NUCLEAR INSTRUMENTS AND METHODS C 

C IN PHYSICS RESEARCH, B5, 74, 1984. C 

C C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INTEGER PCSP 
CHARACTERS PCSZ,PCH 
CHARACTER* 80 XCH , YCH , TCH 
CHARACTER*7 NAME, NAME E 

DIMENSION I A (2048) , IB (2 04 8} , DOPE (12) ,XP (74) ,YP(74) 

DOUBLE PRECISION A(2048) ,B(2048) ,C(2048) ,FR(2048) ,HC(2048) 

DOUBLE PRECISION D(2048) , E (2048) ,G (2) ,S (2048) , SI (2048) 

DOUBLE PRECISION Fl (2048) , PI, V,W,W1 

COMPLEX*16 R(4096),H(4096) ,F(4096) ,T,A1,H1(4096) 

COMPLEX*16 A2,R1(4096) ,Tl 

PI-3.141592653589793D0 

PRINT*, "THIS PROGRAM DEALS WITH THE CONVOLUTION * 
PRINT*, 'H(X)=»F(X) *R(X) , WHERE R IS THE* 
PRINT*, 'RESOLUTION FUNCTION, H IS THE * 
PRINT*, 'BROADENED AND OBSERVED FUNCTION, AND* 
PRINT*, 'F IS THE TRUE AND DESIRED FUNCTION.' 
PRINT*,' ' 

C THE RESOLUTION FUNCTION IS USUALLY OBTAINED FROM 

C THE Sr(85) GAMMA RAY (514 KeV) . 

C 

2 PRINT*, 'THE PROGRAM REQUIRES N DATA POINTS' 
PRINT*, 'WHERE N=2**M, AND M IS AK INTEGER.' 
PRINT* , ' ' 

C THESE DATA POINTS WILL BE H(X) AND R(X) AS COLLECTED 
C BY THE MCA. 

C 

PRINT*, 'ENTER THE VALUE OF H TO BE USED IN THIS' 

PRINT*, 'PROGRAM. (M=ll YIELDS N=2048) » 

PRINT*, ' ' 

READ(*,'(I2) ') M 

N=2**M 

ND2=N/2 

ND4=N/4 

PRINT* , ' ENTER THE RESOLUTION FUNCTION FILENAME' 
PRINT*, ' ' 
READ (*,' (A7) ') NAME 
OPEN (90,FILE=;^AME,PAD= , yES') 
READ (90,1) (IA(I), 1=1,2048) 
1 FORMAT^(5X,5l7) . 
CLOSE (90) 
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PRINT* ' ■ ENTBR • THE "^RIMENTAL FUNCTION FILENAME ' 

READ(*J'(A7) •) NAME 

OPEN (90 , PILE=NAME,PAD s ' YES ' ) 

READ (90,1) (IB(X), 1=1,2048) 

CLOSE (90) 

NIN-2048 

C THE SOBROOTINB DELSKIP DELETES THE LINE 
C SKIP IH THE RAW DATA SUPPLIED BY THE MCA. 

CALL DELSKIP (IA,A,NIN) 

CALL DELSKIP (IB,B,NIN) 
C THE FOLLOWING SOBRODTINE CALCULATES THE CENTRAL 
C CHANNEL LOCATION BASED ON THE CENTER OP THE 
C 11 CHANNEL REGION WHICH CONTAINS THE MAXIMUM 
C TOTAL NUMBER OP COUNTS 

CALL CENTER (A,IMAXA) 

CALL CENTER (B,IMAXB) 

C 

C THIS NEXT LARGE SECTION DEALS WITH THE 
C BACKGROUND CORRECTION. 



C 



PRINT*,* DO YOU WISH TO BACKGROUND CORRECT? » 

READ (*,'(A1) •) NAME 

IP (NAME .EQ. *Y * ) GO TO 3 

CALL NOBACK (N,B,IMAXB) 

GO TO 200 

3 PRINT*, 'DO YOU WISH TO BACKGROUND CORRECT USING* 

PRINT*, 'THE JACKMAN METHOD OR THE SPECTRAL FILTERING* 

PRINT* , * TECHNIQUE? (J OR S) • 

READ (*,* (Al) *) NAME 

IF (NAME .EQ. *J*) GO TO 5 

PRINT*, 'THE NEXT SEVERAL QUESTIONS REFER TO THE* 
PRINT*, 'RESOLUTION FUNCTION. • 
CALL BACK (A,IMAXA,N) 

PRINT*, 'THE NEXT SEVERAL QUESTIONS REFER TO THE' 
PRINT*, 'EXPERIMENTAL FUNCTION.' 
CALL BACK (6, IMAXB, N) 
GO TO 377 

5 PRINT*, 'ENTER THE RESOLUTION OF THE MCA IN' 
PRINT*, 'eV/CHANNEL. (USE A DECIMAL POINT) • 
PRINT*, • (SUGGESTION: IP THE SYSTEM IS SET' 
PRINT*,' IN ITS NORMAL MODE USE 64. )' 
PRINT*,' • 
READ (*,'(F6.2) ') EV 
X-U000.D0/EV) +0.5D0 
IK*IFIX (X) 
IKLA= ( IMAXA-20) /IK 
ILA=514-IKLA 

IRHA=(1024-IMAXA)/IK 

IHA=514+IKHA 

IKLB= ( IMAXB-20) /IK 



CALL PIT (Nl,D,E,G) 

BAH=G(1) 

XMAH=G(2) 



PRINT* , 
PRINT* , 
PRINT* , 
PRINT*, 
READ (* 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
READ (* 



TO WHICH ENERGIES WOULD YOO LIKE TO* 
EXTRAPOLATE THE LINEAR BACKGROUND 1 
CORRECTION? (SUGGESTION: 506 521) ' 



-*) IL,IH 

'ON THE HIGH ENERGY SIDE OF THE PEAK,* 
•WHAT IS THE HIGHEST ENERGY YOO WISH* 
'TO HAVE CONSIDERED FOR DECONVOLDTION? * 
'SUGGESTION: 521' 

*) I END 
IL-IMA7A- (514-IL) *IK 
IH=IMAXVI-(IH-514)*IK ' 
IEND=IMAXA+(IEND-514) *IK 
IEA=IEND 
DO 55 I=IH,IEND 

A ( I) «A( I) - (XMAH*I+BBH) 
IF (A(I) .LT. 0.) A(I)=0.D0 
55 CONTINUE 

XH=XMAH*IB+BAH 

XL=XMAL*IL+BAL 

XM« (XL-XH) / (FLOAT(IL-IH) ) 

XB*XH- ( (XL-XH) / ( FLOAT ( IL- IH) ) ) * IH 

DO 60 I=IMAXA,IB 

A ( I) =A ( I) - (XM*I+XB) 
IF (A(I) .LT. 0.) A(I)=0.D0 
60 CONTINUE 
GO TO 378 

377 CONTINUE 
IEA-IMAXA+109 
IEB-IMAXB+453 

378 CONTINUE 

C 

C THIS IS THE END OF THE BACKGROUND CORRECTION SEQUENCE 
C NEXT WE NEED TO REFLECT THE HIGH ENERGY HALF OF THE 
C BACKGROUND CORRECTED PEAK ABOUT THE CENTRAL POINT 
C TO FORM A SYMMETRIC PEAK. IN ADDITION A ZEROFILL 
C TO 2048 DATA POINTS IS REQUIRED 
C 

CALL ZEROFILL (A, IMAXA, IEA,N) 

CALL ZEROFILL (B, IMAXB, IEB,N) 

M1=IEA-IMAXA+10 

CALL TAPER (A,Ml,N) 

M2=IEB-IMAXB+20 

CALL TAPER (B,M2,N) 

M3-K1 

CALL SMOOTH (A,M3 r N) 
M4=M2 



CALL SMOOTH (B,M4,N) 
C THE NEXT TWO 00 LOOPS NORMALIZE THE AREA 
C OF THE RESOLUTION AND BROADENED SPECTRA. 
V=0.D0 
W-O.DO 
DO 62 1=1, N 
V-V+A(I) 
W=W+B(I) 

62 CONTINUE 
DO 63 I»1,N 

A(I)=A(I)/V 
B(I)«B(I)/W 

63 CONTINUE 
DO 68 I«1,N 

S(I)=B(I) 
68 CONTINUE 

DO 72 1=1, ND2 

J-ND2+I 

V=S(I) 

S(I)=S(J) 

S(J)=V 
72 CONTINUE 

PRINT* , ' DO YOU WISH TO SAVE THE BACKGROUND ' 
PRINT*, 'CORRECTED SPECTRUM? • 
READ (*,'(A1) ') NAME 
IP (NAME .EQ. 'N') GO TO 76 
PRINT*, 'ENTER THE FILENAME. • 
READ (*,' (A7) •) NAME 
OPEN (9 0, FILE-NAME, PAD= 'YES') 
WRITE (90,74) (S(I),I=1,N) 
74 FORMAT (5X,5D22.10) 
CLOSE (90) 
CALL PARAMS (N,S) 

76 CONTINUE 
DO 77 1=1, N 

R(I)=DCMPLX(A(I) ,0.D0) 
H(I)=DCMPLX(B(I) ,0.D0) 

77 CONTINUE 

CALL REORDER (M,H) 
CALL FFT (M f H) 
CALL REORDER (M,R) 
CALL FFT (M,R) 
CALL DECON (M,F,R,H) 

DOPE (1) =150. 
DOPE (2) =150. 
DOPE (3) =950. 
DOPE (4) =650. 
DOPE (5) =0. 
DOPE(6)=0. 
PCSZ='NS' 
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PCSP-1 

c 

C THIS NEXT SMALL SECTION ALLOWS ONE TO VIEW THE POWER 
C SPECTRUM IF SO DESIRED. 

PRINT*, 'WOULD YOU LIKE TO SEE AND PLOT THE POWER SPECTRUM 1 

PRINT*' •° P EXPERIMENTAL DATA AND THE DECONVOLVED DATA.' 

READ (*,'(A1) ') NAME 

IP (NAME . EQ. 'N') GO TO 179 

DO 78 I-1,ND2 

D(I)*H(I)*(CONJG(H(I))) 
78 CONTINUE 

pSnT* , ;' ENT ? R F1LENAME POR TOE ^-S. OP THE EX. DATA. ' 

READ (*,' (A7) •) NAME 
OPEN ( 90, FILE* NAME, PAD=' YES') 
WRITE (90,901) (I,D(I) ,I=1,ND2) 
901 FORMAT (5 (I4,D22 .12) ) 
CLOSE (90) 
CALL REORDER (M,F) 
CALL FFT (M, F) 
T=0.DO 

DO 277 1=1, N 
T=T+F(I) 

277 CONTINUE 

DO 278 1=1, N 
F(I)=F(I)/T 

278 CONTINUE 

CALL REORDER (M,F) 
CALL FFT (M,F) 
DO 178 1=1, ND2 

E(I)=P(I)*(CONJG(F(I))) 
178 CONTINUE 

PRINT* '! ENT f a THE FILENAME POR THE P.S. OF THE DECON. DATA.' 

READ (*,'(A7)«) NAME 

OPEN (90 , FILE= NAME, PAD= 'YES' ) 

WRITE (90,901) (I,E(I) ,I=1,ND2) 

CLOSE (90) 

NPTS»74 

DO 910 1=1,10 

XP (I) =ALOG10 (FLOAT(I) ) 

YP(I)=E(I) 

YP(I)=ALOG10 (YP(I) ) 
910 CONTINUE 

DO 912 1=15,100,5 

J=8+(I/5) 

XP (J) =ALOGl 0 ( FLOAT ( I) ) 
YP(J)=E(I) 
YP(J)=ALOG10 (YP(J) ) 
912 CONTINUE 
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00 914 1-120,1000,20 
J»23+(I/20) 

XP(J)=ALOG10(PLOAT(I)) 
YP(J)«E(I) 
YP(J)=ALOG10(YP( J) ) 
914 CONTINUE 

XP (74) -ALOG10 (1020. ) 
YP(74)«- 10. 

XCH-'LOG OP THE FODRIER COEFFICIENT NO. • 
MX- 3 4 

YCH-'LOG OP THE POWER* 
MY-16 

PRINT*, •ENTER THE TITLE OF THIS PLOT. • 

PRINT*, * • 

READ (*, • (A80) ') TCH 

PRINT*, 'ENTER THE NO. OF CHARACTERS IN THE TITLE • 
READ (*,'(I2)') NT 
PCH-'l ' 

CALL CKSTOP ( 'TK15 * ,9600) 
CALL CKER 

CALL CKGRAF (DOPE,NPTS,XP, YP,XCH,MX, YCH,MY,TCH,NT,PCSZ 
C,PCSP,PCH) 
DO 930 1-1,10 
YP(I)=D(I) 
YP ( I ) -ALOG10 ( YP ( I) ) 
930 CONTINUE 

DO 932 1-15,100,5 
J-8+CX/5) 
YP(J)=D(I) 
YP ( J ) — AL061 0 ( YP ( J ) ) 
932 CONTINUE 

DO 934 1-120,1020,20 
J=23+(l/20) 
YP(J)-D(I) 
YP { J) — ALOG1 0 ( YP ( J) ) 
934 CONTINUE 
PCH- * 3 • 

CALL CKPLT (DOPE, NPTS, XP, YP.PCSZ ,PCSP, PCH) 
CALL CKMA (0.,0.,1) 
READ (*,*) 
179 CONTINUE 

C 

C AT THIS STAGE THE RESOLUTION AND EXPERIMENTAL 
C SPECTRA HAVE BEEN TRANSFORMED INTO FODRIER SPACE. 
C THIS NEXT SMALL SECTION DETERMINES THE NUMBER OF 
C POSITIVE DECREASING FOURIER COEFFICIENTS IN THE 
C DECONVOLUTED RESULT. THIS INFORMATION CAN THEN 
C BE OSED IN THE FILTERING ROUTINES IF SO DESIRED. 



C 



DO 80 1=2, N 
J-I-l 
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D(I)=REAL(F(I)) 

D (J) -REAL (P (J)) 

IP (D(I) .LT. 0.0) GO TO 81 

IP (D(I) .GT. D{J)) GO TO 81 

80 CONTINUE 

81 IPD-I-1 

WRITE (*,82) IPD 

82 C2?2iJ J5X, 'THE NO. OP POSITIVE DECREASING CO EPS. IS ' 15) 
* RINT f 

PRINT* , 'DO TOO WISH TO OSE A COSINE TAPER FILTER' 
PRINT*, 'TO MODIFY THE PODRIER COEFFICIENTS? • 
READ (*, ' (Al) ') NAME 
IF (NAME .EQ. *N') GO TO 189 

PRINT* , ' ENTER THE NUMBER OP THE FIRST COEFFICIENT* 
PRINT*, 'TO BE SET EQUAL TO ZERO. • 
READ (*,*(I3) ■) iz 

PRj^^'ENTER THE^NOMBER OP COEFFICIENTS TO BE TAPERED. • 

IS-IZ-IT 
XS=REAL(IS) 
XZ-REAL(IZ) 
XY=1.0/ (XS-XZ) 
DO 181 I=IS,IZ 
X=REAL(I) 

F(I)=F(I)*0.5*(1-COS(PI*XY*(X-XZ))) 
181 CONTINUE 

J=IZ+1 
K»(N/2)+l 
DO 183 I«J,K 

P(I)»DCMPLX (O.DO,0.DO) 
183 CONTINUE 
R— K- 1 

DO 185 1=2 ,K 
J-N+2-I 
F(J) =F(I) 
185 CONTINUE 

189 PRINT* , ' DO YOU WISH TO USE A BOXCAR FILTER TO' 
PRINT*, 'MODIFY THE FOURIER COEFFICIENTS IN' 
PRINT* , ' ORDER TO REMOVE TEE HIGH FREQ. NOISE? • 
READ (*,'(A1)') NAME 
IF (NAME .EQ. 'N') GO TO 84 

PRINT*, "ENTER THE NUMBER OF COEFFICIENTS TO BE USED* ' 

READ (*,*) J 

J=J+1 

K=N+2-J 

DO 83 I=J,K 

F(I)=DCMPLX(0.DO,0.D0) 

83 CONTINUE 

84 CONTINUE 

PRINT*, 'DO YOU WANT TO SEE THE FOURIER COEFFICIENTS^' 
READ (*,'(A1)') NAME 



IP (NAME .EQ. *N* ) GO TO 85 
PRINT* , • ENTER THE FILENAME FOR H(X). ■ 
PRINT* , • ' 
READ (*, ' (A7) ') NAME 
OPEN (90, PILE* NAME , PAD= * YES ' ) 
WRITE (90,911) (I,H(I), 1=1, N) 
911 PORMAT (2(14,2022.10)) 
CLOSE (90) 

OPEN (90,FILE='FCSF' ,PAD='YES') 
WRITE (90,911) (I,F(I), 1=1, N) 
CLOSE (90) 

85 CONTINUE 

CALL REORDER (M,P) 

CALL PPT (M,F) 

CALL REORDER (M,B) 

CALL PPT (M,B) 
C THIS NEXT SECTION ALLOWS ONE TO USE TEE ITERATION 
C ROUTINE TO IMPROVE P (X) . 

PRINT* , ' WOULD YOU LIKE TO USE THE ITERATION' 

PRINT*, 'SUBROUTINE TO TRY TO IMPROVE THE' 

PRINT*, 'ACCURACY OP THE DECONVOLUTED SPECTRUM? 

READ (*, ' (Al) ') NAME 

IP (NAME .EQ. 'N') GO TO 86 

CALL ITER (M,R,F,B) 

86 CONTINUE 

C THIS NEXT DO LOOP SHIFTS THE PHASE OF THE 
C EXPERIMENTAL FUNCTION B, AND TBE DECONVQLDED 
C FUNCTION F, SO THAT THEIR MAXIMA ARE LOCATED 
C IN THE CENTRAL CHANNELS. 
ND2=N/2 
DO 88 1=1, ND 2 
J=ND2+I 
T=F(I) 
T1=B(I) 
F(I)=F(J) 
B(I)=B(J) 
F(J)=T 
B(J)=Tl 
88 CONTINUE 

DO 98 1=1, N 

Fl(I)=REAL(F(I) ) 

98 CONTINUE 
W=0.D0 

DO 99 1=1, N 
W=W+P1(I) 

99 CONTINUE 

DO 100 1=1, N 

F1(I)=F1(I)/W 
100 CONTINUE 

PRINT*,'DO YOU WANT A PLOT OF h(X) AND f(X)? ' 
READ (*, ' (Al) ') NAME 
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IP (NAME .EQ. 'N 1 ) GO TO 105 
NPTS-60 
W«P1(1025) 
DO 102 1=1,60 
J-1024+I 

XP(I)=.064*(FLOAT(I)) 
YP(I)=F1(J)/W 

102 CONTINUE 
XCH-'KeV* 
MX-3 

YCB«*h(X) AND f (X) • 
MY=13 

PRINT* , * ENTER THE TITLE OP THIS PLOT. * 

PRINT*, • • 

READ (*,*(A80)») TCH 

PRINT*, •ENTER THE NO. OP CHARACTERS IN THE TITLE. » 
READ (*,*(I2)») NT 
PCH« ' 1 1 

CALL CKSTOP ( 1 TK15 * ,9600) 
CALL CKER 

CALL CKGRAP (DOPE,NPTS,XP,yp, XCH,MX, YCH,MY,TCH,NT,PCSZ 
C,PCSP,PCE) 
DO 103 1=1,60 

J=1024+I 

YP(I)=S(J)/W 

103 CONTINUE 
PCH='3 ' 

CALL CKPLT (DOPE,NPTS,XP,YP,PCS2,PCSP,PCH) 

CALL CKMA (0.,0.,1) 

READ (*,*) 
105 PRINT*, , DO YOO WISH TO SAVE THE ' 

PRINT*, ' DECONVOLUTED SPECTRDM? ' 

READ (*,'(Al) f ) NAME 

IP (NAME .EQ. 'N') GO TO 101 

PRINT*, •ENTER THE FILENAME. * 

READ (*,MA7)») NAME 

OPEN (90, FILE=NAME , PAD= ' YES ' ) 

WRITE (90,74) (Pl(I) ,I=1,N) 

CLOSE (90) 
101 CONTINUE 

PRINT*,* DO YOU WISH TO CHECK YOUR RESULT BY' 

PRINT* , * CONVOLUTING f(X) WITH r(X) AND COMPARING* 

PRINT*, •THE RESULT WITH h(X)?* 

READ (*, * (Al) •) NAME 

IP (NAME .EQ. *N' ) GO TO 200 

CALL REORDER (M,F) 

CALL FFT (M,F) 

CALL CONV (M,F,R,H1) 

CALL REORDER (M,H1) 

CALL FPT (M,H1) 

DO 106 1=1, N 
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BC(I)«REAL(Hl(I)) 

106 CONTINUE 
V*0 . DO 

DO 107 1=1, N 
V=V+HC(I) 

107 CONTINUE 

DO 108 1*1, N 

HC(I)=HC(I)/V 

108 CONTINUE 

PRINT* , 1 ENTER TEE FILENAME POR THE RESULT ' 

PRINT*,* OP THIS COMPARISON.' 

PRINT*,* * 

READ (*,'(A7)') NAME 

OPEN {90 , PILE=NAME, PAD= * YES ' ) 

WRITE (90,110) 

110 FORMAT (3X, , I , ,14X, , h(X)RAW , ,14X,»h(X)»(F*R) »,14X,'f (X) ') 
WRITE (90,120) (I, S(I) ,HC(I) ,F1(I) ,1=1025,1300) 

120 FORMAT (1X,I4,2X,3D22.10) 
CLOSE (90) 

PRINT*, 'DO TOU WANT A PLOT OF h(X) AND THE ESTIMATE' 
PRINT*, ' OF h(X) OBTAINED FROM f (X)*r(X)? 1 
READ (*, * (Al) ') NAME 
IF (NAME -EQ. *N ( ) GO TO 200 
KPTS=60 
W=HC (1025) 
DO 130 1=1,60 
J=1024+I 

XP(I)=.064*(FLOAT(I)) 

TP(I)«HC(J)/W 
130 CONTINUE 
XCB='KeV» 
MX=3 

YCE='h(X) AND h* (X) 1 
MY-14 

PRINT*, 'ENTER THE TITLE OF THIS PLOT.' 

PRINT* , ' ' 

READ (*,' (A80) ') TCH 

PRINT* , * ENTER THE NO. OF CHARACTERS IN THE TITLE. • 
READ (*. ' (12) ') NT 
PCH-'l ' 

CALL CKSTUP ( 'TK15 ' ,9600i 
CALL CKER 

CALL CKGRAF (DOPE,NPTS,XP, YP,XCH,MX, YCB, MY,TCB,NT,PCSZ 
C,PCSP,PCH) 
DO 140 1=1,60 
J=1024+I 
YP(I)=S(J)/W 
140 CONTINUE 
PCH='3 ' 

CALL CKPLT (DOPE,NPTS,XP, YP, PCSZ ,PCSP, PCH) 
CALL CKMA (0. ,0. ,1) 



READ (*,*) 
200 CONTINUE 

CALL P ARAMS (N,P1) 

STOP 

END 

cccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE DELSKIP (IC,C,NIN) 
DIMENSION IC(2048) 

DOUBLE PRECISION C(1024) 
J=0 

K-0 

DO 10 I«1,NIN 
J=J+1 

IC(I)«IC(K) 
IP (J .EQ. 5) THEN 
J=0 
K»K+5 
END IP 
10 CONTINUE 
NID2=NIN/2 
DO 20 I«=1,NID2 

C(I)«PLOAT(IC(I)) 
20 CONTINUE 
RETURN 
END 

C 

SUBROUTINE CENTER (A, I MAX) 

DOUBLE PRECISION A (102 4) ,X,T(21) ,Y 

X=0 .DO 

DO 10 1=300,800 

IP (A (I) .GT. X) THEN 

X-A(I) 

J-I-10 

ENDIF 
10 CONTINUE 
K-J+21 
Y-0.D0 
DO 20 I=J,K 

L=I-5 

M=I+5 

T(I)~O.D0 

DO 15 N»L,M 

T(I)«T(I)+A(N) 
15 CONTINUE 

IP (T(I) .GT. Y) THEN 

Y=T(I) 

I MAX* I 

ENDIF 
20 CONTINUE 
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RETURN 
END 

C 

SUBROUTINE FIT (N,X,Y,C) 
C THIS SUBROUTINE IS A SIMPLIFICATION OF 
C A LARGER SUBROUTINE WBICH FITS A CURVE USING 
C UP TO A NINTH ORDER POLYNOMIAL. FOR A 
C REFERENCE SEE JPS DISSERTATION. 

DOUBLE PRECISION X(300) ,Y(300) ,C(2) ,A(3,3) f XN(300) 
DO 10 I«1,N 
10 XN(I)=1.D0 
DO 30 1-1,2 

A(I,1)«=0.D0 
A(I,3)»O.D0 
DO 20 J=1,N 

A(I,1)=A(I,1)+XN(J> 
A(I,3)-A(I,3)+Y(J)*XN(J) 
20 XN(J)=XN(J)*X(J) 
30 CONTINUE 

A(2,2)=0.D0 
DO 40 J=1,N 

A(2,2)=A(2,2)+XN(J) 
40 XN(J)=XN(J)*X(J) 
A(1,2)=A(2,1) 
DO 50 1-1,2 
SUM=0.D0 

IF (2 .GT. I) GO TO 45 
SUM=SUM+A (1,1) *A(1,2) 
A(I,2)=A(I,2) -SUM 
GO TO 50 

45 IF (ABS(A(I,I)) .LT. l.E-10) GO TO 99 

A(I,2)«(A(I,2)-SUM)/A(I,I) 
50 CONTINUE 

DO 90 J=l r 2 
90 C(J)=A(J,3) 

C(1)«C(1)/A(1,1) 

SUM=A(2,1)*C(1) 

C(2)«(C(2)-SUM)/A(2,2) 

SUM=A(1,2)*C(2) 

C(l)=C(l)-SOM 

RETURN 

99 PRINT*, 'TROUBLE IN SUBROUTINE FIT' 
RETURN 
END 

C 

SUBROUTINE ZEROFILL (C,IM,IE,N) 
DOUBLE PRECISION C (2048) ,D (2048) , T 
M=IE-IM 
MP1=M+1 
DO 10 1*1, N 
D(I)=0.D0 



10 CONTINUE 

DO 20 1=1, MP 1 
J-IM+I-1 
D<I)«C(J) 
20 CONTINUE 

DO 30 I«1,N 
C(I)=D(I) 
30 CONTINUE 
RETURN 
END 

SUBROUTINE TAPER (C,M,N) 
DOUBLE PRECISION C(2048),T 
T-l. DO/1 000. DO 
DO 90 I«1,M 

IF (C(I) .NE. 0.DO) GO TO 80 

J=I-1 

K*I+1 

C(I) = (C(J)+C(K))/2.D0 

IP (C(I) .LT. T) GO TO 95 

80 CONTINUE 

90 CONTINUE 
GO TO 97 

95 IK»I 

DO 96 J«IK,M 
K=J-1 

C(J)=0.9 5D0*C(K) 

96 CONTINUE 

97 CONTINUE 
ND2=N/2.D0 

DO 100 1=2, ND2 

J=ND2+I 

K«ND2+2-I 

C(J) =C (K) 
100 CONTINUE 
RETURN 
END 

SUBROUTINE SMOOTH (C f M,N) 
DOUBLE PRECISION C(2048) 
DO 20 1=2, M 
J=I-1 

IP (C(J) .GT. C(I)) GO TO 10 
K=I+1 

C(I)=(C(J)+C(K))/2.D0 
5 IF (C(J) .GT. C(I)) GO TO 10 
C(I)=0.999*C(I) 
GO TO 5 
10 CONTINUE 
20 CONTINUE 
ND2=N/2 
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DO 30 1-2, ND2 

J-ND2+I 

K-ND2+2-X 

C(J)-C(K) 
30 CONTINUE 
RETURN 
END 

C 

SUBROUTINE REORDER (M,C) 
C IN ORDER TO USE THE FOLLOWING FFT ONE MUST 
C FIRST REORDER THE DATA. SEE JPS DISSERTATION 
C POR DETAILS. 

COMPLEX*16 C(4096) ,T 

N«2**M 

MV2«=N/2 

NM1=N-1 

J=l 

DO 7 1=1, NM1 

IF (I .GE. J) GO TO 5 

T=C(J) 

C(J)-C{I) 

C(I)=T 

5 K=MV2 

6 IF (K .GE. J) GO TO 7 
J=J-K 

K=K/2 
GO TO 6 

7 J=J+K 
RETURN 
END 



SUBROUTINE FFT (K, A) 
C FOR DETAILS SEE JPS DISSERTATION 
DOUBLE PRECISION PI 
COMFLEX*16 A(4096) ,U,W,T 
PI»3.141592653589793D0 
N«2**M 
DO 30 L=1.M 

LE=2**L 

LEI=LE/2 

U= (l.D0,0.D0) 

W=DCMPLX ( DCOS ( PI/LEI ) , -DS IN ( PI/LE I ) ) 
DO 20 J=1,LEI 

DO 10 I=J,N,LE 
IP=I+LEI 
T=A(IP) *D 
A(IP)=A{I)-T 
A(I)=A(I)+T 
10 CONTINUE 
U=U*W 



20 CONTINUE 
30 CONTINUE 

RETURN 

END 

C 

SUBROUTINE CONV (M,F,R,H) 
DOUBLE PRECISION X 

COMPLEX*16 B(4096) ,F(4096) ,R(4096) 

N=2**M 

Wl«FLOAT(N) 

DO 10 1=1, N 

B(I)-P(I)*R(I) 
' 10 CONTINUE 

DO 20 I«1,N 

H(I)«H(I)/Wl 
20 CONTINUE 
RETURN 
END 

C 

SUBROUTINE DECON (M,F,R,H) 
DOUBLE PRECISION X 

COMPLEX*16 F(4096) ,R<4096) ,H(4096) 
N=2**M 

C IN FOURIER SPACE DECONVOLUTION IS SIMPLT A 
C COMPLEX DIVISION. 
DO 10 1*1, N 

F(I)=E(I)/R(I) 
10 CONTINUE 
X* FLOAT (N) 
DO 30 1=1, N 

F(I)=F(I)/X 
30 CONTINUE 
RETURN 
END 

C 

SUBROUTINE ITER (M,R,F,H) 
DIMENSION A(2048) ,B(2048) 

COMPLEX*16 R(4096) ,F(4096) ,H(4096) ,H1(4096) ,T 

CHARACTER* 7 NAME 

PRINT*, 'ENTER THE VALUE FOR THE PARAMETER WHICH ' 

PRINT*, •CONTROLS THE ALLOWABLE ERROR BETWEEN THE 

PRINT*, 'OBSERVED SPECTRUM AND F*R. • 

PRINT* , ' EPS=THE SQRT OF THE SUM OF THE SQRS OF' 

PRINT*, 'THE DIFFERENCES BETWEEN THE H(I)S AND' 

PRINT*, 'THE F(I)*R(I)S. 

PRINT* , * (A REAL NUM3ER PLEASE.) ' 

READ (*, ' (F13.7) ') EPS 

PRINT* , ' ENTER THE MAXIMUM NUMBER OF ITERATIONS' 
PRINT* , ' YOU WISH TO PERMIT. ' 



irtixNi *' , ENTER THE MAXIMUM N 
PRINT* , ' YOU WISH TO PERMIT. 
PRINT*, ' (AN INTEGER PLEASE) 
READ (*, ' (13) ') I END 
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N«2**M 

ND2-N/2 

ND4-N/4 

DO 30 I-1,IEND 

CALL REORDER (M,F) 

CALL PPT (M,P) 

CALL CONV (M,F,R,B1) 

CALL REORDER (M,F) 

CALL FPT (M,F) 

CALL REORDER (M,B1) 

CALL PPT <M,B1) 

RSQ»0.0 

DO 2 IA=1,N 

A(IA)=REAL(B(IA) ) 

B(IA)-R£AL(Bl(IA)) 
2 CONTINUE 
V=0.0 
W=0.0 

DO 4 IAA=1 , N 
V-V+A(IAA) 
W-W+B(IAA) 
4 CONTINUE 

DO 6 IAB«1,N 

A(IAB)=A(IAB)/V 
B(IAB)=B(IAB)/W 
6 CONTINUE 

DO 10 IAC=1,ND4 

RSQ=RSQ+ ( A ( IAC) -B ( IAC) ) * *2 

10 CONTINUE 
RSQ=RSQ**0.5 
WRITE (*,11) RSQ 

11 FORMAT <5X,»RSQ« »,E13.6) 

IP (RSQ -LT. EPS) GO TO 33 

P(1)=P(1)*(A(1)+2*A(2)+2*A(3) )/ (B (1) +2*B (2) +2*B(3) ) 

F(2)~F<2)*(A(l)+2*A(2)+A(3)-t-A{4))/(B(l)+2*B(2)+B(3)+B(4)) 
DO 20 IB-3,ND4 
J=IB-2 
K=IB+2 
TB-0.0 
TH1-0.0 
DO 15 II=J,K 
TH«TB+A(II) 
TH1=TH1+B(II) 
15 CONTINUE 

IF (TBI . EQ. 0.0) GO TO 20 
F ( IB) «=F { IB) *TB/TBl 
20 CONTINUE 
DO 25 IC-2,ND2 
J=ND2+IC 
K=ND2+2-IC 
F(J)=F(K) 



25 CONTINUE 
30 CONTINUE 
33 CONTINUE 
I-I-l 

WRITE (*,35) I 

(5X ' ,THE N °- ° P VENATIONS IS M3, 

END 

C 

SUBROUTINE P ARAMS (N,A) 

3SE££g?E£ A(2048 >' x <").M(io,,v, W 

K=1024 

DO 20 1=1,12 
X(I)«0.D0 
DO 10 J=l,5 
K=K+1 

in ~ *<I)«X(I)+A(K) 

10 CONTINUE 
20 CONTINUE 
V=0.D0 

DO 30 1=1025,1084 
V=V+A(I) 
30 CONTINUE 

DO 40 1=1,12 

40 CONTINUE 

DO 60 1=2,10,2 
M(I)=0.D0 
DO 50 J=1026,1084 
V=PLOAT(J-1025) 
W=FLOAT(I) 
V=V**W 

_ M(D-M(I)+A(J)*V 
50 CONTINUE 

M(I)=M(I)/59.D0 
60 CONTINUE 

M(2)=OSQRT(M(2) ) 

V*M(2) 

DO 65 1=4,10,2 

„ _ M(D=M(I)/(V**I) 
65 CONTINUE 

PRINT*, 'ENTER THE FILENAME FOR trtc orm «„, 
PRINT*, -PARAMETERS. « * TOIS SET OP 

READ (*,»(A7) •) NAME 

OPEN (90,FILE=NAME,PAD=* YES' ) 
WRITE (90,70) ' 

70 FORMAT (5X, 'EVEN CENTRAL MOMENTS') 

*n W J£ TE (90 ' 80) <^M(I,ri=2,8,2 f } 
80 FORMAT (10X, 'M( ' ,11, • ) =',D1710) 
WRITE (90,100) (M(10)) fU1, ' 1 °' 



100 FORMAT (9X,'M(10) «',D17.10) 

WRITE (90,110) 
110 FORMAT (5X,'5 CHANNEL GROUPS (CENTER FIRST! 

WRITE (90,120) (I,X(I),I-1,9) 
120 FORMAT (10X, »X( • ,11, ') =*,D17.10) 

WRITE (90,130) (I, X(I) ,1=10,12) 
130 PORMAT (9X, 'X(',I2,') = ',D17.10) 

CLOSE (90) 

RETURN 

END 

SUBROUTINE BACK (A,IMAX,N) 
DOUBLE PRECISION A (1024) ,C (2048) , V 
COMPLEX *1 6 B(4 096) 
CHARACTER* 7 NAME 
ND2=N/2 
ND4=N/4 
DO 10 1-1, N 
C(I)=0.D0 
10 CONTINUE 

N1=ND2-IMAX 

IF (Nl .6T. ND4) N1=ND4 
DO 20 1=1, Nl 
J=IMAX+I-1 
C(I)=A(J) 
20 CONTINUE 
N2=IMAX-1 
N3=N-N2 
DO 30 J=1,N2 
I=N3+J 
C(I)-A(J) 
30 CONTINUE 

DO 40 1=1, N 

B(I)=DCMPLX(C(I) ,0.D0) 

40 CONTINUE 
M=ll 

CALL REORDER (M,B) 
CALL FFT (M,B) 

PRINT*, 'WOULD YOU LIKE TO SEE THE FIRST FEW* 

PRINT*, 'FOURIER COEFFICIENTS? • 

READ (*, • (Al) •) NAME 

IF (NAME .EQ. 'N') GO TO 43 

PRINT* , • HOW MANY? • 

READ (*,*) IFC 

WRITE (*,41) (I,B(I),I«1,IFC) 

41 FORMAT (I4,2D22.10) 
43 CONTINUE 

PRINT*, 'WOULD YOU LIKE TO SEE THE FIRST FEW* 
PRINT*, 'VALUES OF THE POWER SPECTRUM? ' 
READ (*, • (Al) ') NAME 
IF (NAME . EQ. 'N') GO TO 46 



44 



PRINT* , • HOW MANY? • 
READ (*,*) IPS 
DO 44 I«1,IPS 

A(I)*B(I)*CONJG(B(I)) 
- . CONTINUE 

WRITE (*,45) (I,A(I),I=1,IPS) 

45 FORMAT (I4,D22.10) 

46 CONTINUE 

PRINT*,' DO YOU WISB TO DELETE ALL OP THE SINE ' 
PRINT* , ' COEFFICIENTS? ' 
READ (*, ' (Al) •) NAME 
IF (NAME . EQ. 'N') GO TO SI 
DO 50 1=1, ND 2 

B(I)«(B(I)+CONJG(B(I) ) )/2.D0 

50 CONTINUE 

51 PRINT* , • DO YOU WISH TO DELETE SOME OF THE COSINE* 
PRINT*, 'COEFFICIENTS? ' 
READ {*, • (Al) •) NAME 
IF (NAME .EQ. 'N') GO TO 53 
PRINT*, 'HOW MANY? ' 
READ (*,*) IDC 
DO 52 1=1, IDC 

B(I)»(B(I)-CONJG(B(I) ))/2.D0 

52 CONTINUE 
DO 152 1*2, IDC 

J=2050-I 

B(J)=(B(J)-CONJG(B(J) ) )/2.D0 
.. CONTINUE 

53 PRINT* ,' INSTEAD OF DELETING SOME OF THE COSINE' 
PRINT*, 'COEFS. YOU MAY PREFER TO JUST TRIM THEM. " 
PRINT*, 'THIS IS DONE BY SETTING THE FIRST N TERMS 
PRINT* ,' EQUAL TO THE N+l TERM. YOU MAY INVOKE ' 
PRINT*, 'THIS OPTION EVEN IF YOU HAVE PREVIOUSLY' 
PRINT*, 'DELETED THE FIRST N COSINE COEFFICIENTS.' 
PRINT*, 'DO YOU WISH TO DO SO? ' 

READ (*,'(A1)») NAME 
IF (NAME .EQ. »N') GO TO 200 
PRINT*, 'BOW MANY? • 
READ (*,*) IT 
DO 190 1=1, IT 
J=IT+1 

B(I)=DCMPLX(REAL(B(J)) ,0.D0) 
190 CONTINUE 
200 CONTINUE 

PRINT*, 'DO YOU WISH TO DELETE SOME OF THE SINE' 

PRINT*, 'COEFFICIENTS? ' 

READ (*,' (Al) ') NAME 

IF (NAME . EQ. 'N' ) GO TO 55 

PRINT* , * HOW MANY? • 

READ (*,*) IDS 

DO 54 1=1. IDS 
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B(I) = (B(I)+CONJG(B(I)))/2.D0 

54 CONTINUE 

DO 154 1=2, IDS 
J-2050-I 

B(J) = (B(J) +CONJG(B(J)) )/2.D0 
154 CO NT IN DE 

55 CONTINUE 

PRINT*,' DO YOU WISH TO SAVE THE FOURIER COEFS. OF' 

PRINT*, *B (X) PRIOR TO THE TAPER IN THE WINGS? ' 

READ (*,' (Al) ') NAME 

IF (NAME .EQ. 'N') GO TO 57 

PRINT*, 'ENTER THE FILENAME. ' 

PRINT*,' ' 

READ (*,'(A7)') NAME 

OPEN (90, FILE-NAME, PAD* 'YES') 

WRITE (90,56) (I ,B ( I) ,1=1,2048) 

56 FORMAT (2 (14 ,2D22 .10) ) 
CLOSE (90) 

57 CALL REORDER (M,B) 
CALL FFT (M,B) 

DO 60 1=1, ND2 

C(I)=REAL(B(I)) 
60 CONTINUE 
V«O.D0 

DO 70 1=359,484 
V=V+C(I) 
70 CONTINUE 
V-V/126. 

PRINT* , ' DO YOU WANT A D.C. CORRECTION USING A' 
PRINT*, 'LEAST SQUARE FITTING IN THE TAILS? ' 
READ (*, 1 (Al) ') NAME 
IF (NAME .EQ. 'N') V=0.D0 
DO 80 1=1, ND 2 
C(I)=C(I)-V 

IF (C(I) .LT. 0.D0) C(I)=0.D0 
80 CONTINUE 

DO 90 1=1, ND4 
J=(3*ND4)+I 
A(I)=C(J) 
90 CONTINUE 
N4=ND4+1 
DO 100 I=N4,ND2 
J=I-ND4 
A(I)=C(J) 
100 CONTINUE 
IMAX=513 
RETURN 
END 

SUBROUTINE NOB AC K (N,B,IMAX) 
DOUBLE PRECISION B (1024) ,A (2048) , V 



C0MPLEX*16 C(4096) 
CHARACTER* 7 NAME 
M»ll 
ND2=N/2 
ND4=N/4 
DO 10 1=1, N 
A(I)=O.D0 
10 CONTINUE 

N1=ND2-IMAX 

IF (Nl .GT. ND4) N1=ND4 
DO 20 1*1, Nl 
J=IMAX+I-1 
A(I)«B(J) 
20 CONTINUE 
N2-IMAX-1 
N3=N-K2 
DO 30 J=1,N2 
I=N3+J 
A(I)«B(J) 
30 CONTINUE 
V=0.D0 
DO 33 1=1, N 
V=V+A(I) 
33 CONTINUE 

DO 36 1=1, N 

A(I)=A(I)/V 
36 CONTINUE 

OPEN (90 # FILE=' RAWDATA' ,PAD= , yES' ) 
WRITE (90,70) (I,A(I),I«1,N) 
CLOSE (90) 
DO 40 1=1, N 

C(I)=DCMPLX(A(I) ,0.D0) 
40 CONTINUE 

CALL REORDER (M,C) 
CALL FPT (M,C) 
DO 50 1=1, N 

A(I)=CiI)*CONJG(C(I) ) 
50 CONTINUE 

PRINT* , ' ENTER THE FILENAME FOR E (X) » 
PRINT* , ' ' 
READ {*, • (A7) ») NAME 
OPEN ( 90, FILE-NAME, PAD= 'YES' ) 
WRITE (90,60) (I, C(I) ,1=1,2048) 
60 FORMAT (2(I4,2D22.10)) 
CLOSE (90) 

PRINT*, 'ENTER THE PLENA ME FOR THE P.S. 
PRINT* , ' ' 
READ (*, ' (A7) ') NAME 
OPEN (90 , FILE-NAME,PAD= ' YES' ) 
WRITE (90,70) (1,A(I) ,1=1,1024) 
70 FORMAT (5 ( I4,D22 .10) ) 



CLOSE (90) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccc 



APPENDIX 6 



In this Appendix the theoretical relationship between the vacancy 
concentration, C v , and the concentration of dislocation cores, C d , is 
demonstrated, (from CotterHl et aK 1972). The starting equation is 

C v = A exp(-E vf /KT) (6 . 1} 
where C v 1s the vacancy concentration at thermal equilibrium 

E v f is the vacancy formation energy 
and KT has its usual meaning 

Now the concentration of vacancies after the quench, C*. is related to C by 

v v ' 



C v = nC v (G.2) 
where n is the quench efficiency which is assumed to be nearly unity 

After the low temperature anneal the concentration of dislocation 
loops, C t , is given by 

C A " vC v (6.3) 
where v is the inverse of the number of vacancies required to form an 
average loop 

The dislocation concentration, C d , in terms of an atomic fraction is then 

C d " <VD/d 0 (6 .4, 
where 0 is the average loop diameter 
and d 0 is the nearest neighbor distance. 
It can be shown that 
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« - 2V~3 d 2 /itD 2 (6.S) 
Combining Equations (G.3), (6.4), and (G.5) yields 

C d " 2 ^ d o C > (6-6) 



APPENDIX H 



In this Appendix the derivation of the defect specific lineshape 
parameter M R" is recreated from the work of Mantl and Trifthauser (1978). The 
following quantities require definition: 

S,W - are the values of the tt S" and "W M parameters for the sample 
being analyzed. 

S t ,W t - are the values of the »S U and M W" parameters for a sample in 
which all of the positrons are annihilating from the trapped 
state. 

Sb.W b - are the values of the M S W and M W M parameters for a sample in 
which all of the positrons are annihilated from the bulk 
(defect free). 

P t - the probability of a positron annihilation from the trapped 
state. 

P B - the probability of a positron annihilating from the bulk. 
m - the trapping constant. 
C - the trap concentration. 

X B - the annihilation rate from the bulk state. 
Now, it is known that 

S = S B P B + S t P t (H# i) 

and 

P B + Pt = 1 (H.2) 
(143) 
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AlSO P t = yC/(X B + uC) (H.3) 
Substituting Eq. (H.3) Into Eq. (H.l) yields: 

s = h » - v?Sc> ♦ * t ij^i ,„. 4 , 

Similarly 

W = "b * ' A-5nrc> + « t (H.5) 

Rearranging Eq. (H.4) and Eq. (H.5) yields: 

S - S B = (S t - S B )pC/(A f + mC) (h.6) 

W - W B = (W t - W B ) M C/(X f + uC) (H.7) 

Taking the ratio of Eq. (H.6) to Eq. (H.7) yields: 



S 



s t -s. 



t b| 

Clearly the ratio on the far right has a constant value for a given defect in 
a given material . 
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